Covariation of waking cortisol with DHEA and testosterone across the pubertal transition: Effects of threat-related early life stress

#Libraries
library(papaja)
library(tidyverse)
library(lme4)
library(lmerTest)
library(modelr)
library(sjstats)
#library(olsrr)
library(broom)
library(haven)
library(corrr)
library(ggcorrplot)
library(corrplot)
library(cowplot)
#Files 
cd_lf_all_file <- "~/Box/lucy_king_files/ELS/cort_dhea/data/data_final/cort_dhea_all_20190509.csv"
cd_lf_file <- "~/Box/lucy_king_files/ELS/cort_dhea/data/data_final/cort_dhea_final_20190509.csv"
cd_wf_all_file <- "~/Box/lucy_king_files/ELS/cort_dhea/data/data_final/cort_dhea_all_wf_20190509.csv"
cd_wf_file <- "~/Box/lucy_king_files/ELS/cort_dhea/data/data_final/cort_dhea_final_wf_20190509.csv"
els_file <- "~/Box/lucy_king_files/ELS/stress_data/data/els_severity_scores_20190221.csv"
severity_file <- "~/Box/lucy_king_files/ELS/stress_data/data/severity_by_event_20190825.csv"
demo_file <- "~/Box/lucy_king_files/ELS/cort_dhea/data/data_final/demo.csv"
tanner_file <- "~/Box/lucy_king_files/ELS/cort_dhea/data/data_final/tanner_t1t2.csv"
#Themes
theme_set(theme_apa())
#Functions
format_pval <- function(pval, digits = 3){
  pval <- str_replace(round(signif(pval), digits), "0.", ".")
  pval <- if_else(pval == 0, "<.001", paste0("=", pval))
  return(pval)
}
format_lm_df <- function(mod) {
  df <- broom::tidy(mod)
  
  df %>% 
    gather(key, value, estimate:p.value) %>% 
    unite(new, term, key) %>% 
    spread(new, value)
}
set.seed(123)
options(scipen = 999)

cd_wf_all <-
  read_csv(cd_wf_all_file)

Parsed with column specification:
cols(
  .default = col_double(),
  session_type = col_character(),
  timepoint = col_character(),
  scan_date = col_datetime(format = ""),
  date_collection = col_date(format = ""),
  medication = col_character(),
  birth_date = col_datetime(format = ""),
  tanner_date = col_date(format = "")
)
See spec(...) for full column specifications.

cd_wf_clean <-
  read_csv(cd_wf_file) %>% 
  mutate(
    male = as.factor(male),
    medication_t1 = as.factor(medication_t1),
    medication_t2 = as.factor(medication_t2)
  )

Parsed with column specification:
cols(
  .default = col_double(),
  session_type = col_character(),
  timepoint = col_character(),
  scan_date = col_datetime(format = ""),
  date_collection = col_date(format = ""),
  medication = col_character(),
  birth_date = col_datetime(format = ""),
  tanner_date = col_date(format = "")
)
See spec(...) for full column specifications.

cd_all <- 
  read_csv(cd_lf_all_file) %>% 
  mutate(
    timepoint = as.factor(timepoint)
  )

Parsed with column specification:
cols(
  .default = col_double(),
  session_type = col_character(),
  timepoint = col_character(),
  scan_date = col_datetime(format = ""),
  date_collection = col_date(format = ""),
  medication = col_character(),
  birth_date = col_datetime(format = ""),
  tanner_date = col_date(format = "")
)
See spec(...) for full column specifications.

cd_clean <- 
  read_csv(cd_lf_file) %>% 
  mutate(
    timepoint = as.factor(timepoint),
    male = as.factor(male),
    medication_binary = as.factor(medication_binary)
  )

Parsed with column specification:
cols(
  .default = col_double(),
  session_type = col_character(),
  timepoint = col_character(),
  scan_date = col_datetime(format = ""),
  date_collection = col_date(format = ""),
  medication = col_character(),
  birth_date = col_datetime(format = ""),
  tanner_date = col_date(format = "")
)
See spec(...) for full column specifications.

# create person-mean-centered variable for tanner stage and average across t1 and t2 for other time-varying covariates
cd_wf_clean <-
  cd_wf_clean %>% 
  mutate(
    tan_diff = if_else(
      !is.na(tanner_t2) & !is.na(tanner_t1),
      tanner_t2 - tanner_t1, NA_real_
    ),
    dhea_diff = if_else(
      !is.na(dhea_clean_t2) & !is.na(dhea_clean_t1),
      dhea_clean_t2 - dhea_clean_t1, NA_real_
    ),
    test_diff = if_else(
      !is.na(test_clean_t2) & !is.na(test_clean_t1),
      test_clean_t2 - test_clean_t1, NA_real_
    )
  )
cd_clean <- 
  cd_clean %>% 
  left_join(
    cd_wf_clean %>% 
      filter(
        # prevents from merging a difference score for those with hormones at only 1 timepoint (do not want to model their change)
        !is.na(cort_clean_t2) | !is.na(dhea_clean_t2) | !is.na(test_clean_t2),
        !is.na(cort_clean_t1) | !is.na(dhea_clean_t1) | !is.na(test_clean_t1)
      ) %>% 
      select(ELS_ID, tan_diff, dhea_diff, test_diff), 
    by = "ELS_ID"
  ) %>% 
  group_by(ELS_ID) %>% 
  mutate(
     time = case_when(
      !is.na(interval_yr) & timepoint == "T1" ~ 0,
      !is.na(interval_yr) & timepoint == "T2" ~ interval_yr,
      TRUE ~ NA_real_
    ), 
    time_per_cent = if_else(
      !is.na(interval_yr),
      time - mean(time), NA_real_
    ),
    dhea_per_cent = if_else(
      !is.na(dhea_diff),
      dhea_clean - mean(dhea_clean), NA_real_
    ),
    test_per_cent = if_else(
      !is.na(test_diff),
      test_clean - mean(test_clean), NA_real_
    ),
    tan_per_cent = if_else(
      !is.na(tan_diff),
      tanner_av - mean(tanner_av), NA_real_
    ),
    tan_diff_lg = case_when(
      timepoint == "T1" & !is.na(tan_diff) ~ 0,
      timepoint == "T2" ~ tan_diff,
      TRUE ~ NA_real_     
    ),
    age_av_t1t2 = mean(child_age, na.rm = TRUE),
    dhea_av_t1t2 = mean(dhea_clean, na.rm = TRUE),
    test_av_t1t2 = mean(test_clean, na.rm = TRUE),
    tanner_av_t1t2 = mean(tanner_av, na.rm = TRUE),
    timecoll_av_t1t2 = mean(time_collection, na.rm = TRUE),
    bmi_av_t1t2 = mean(bmi, na.rm = TRUE),
    age_av_t1t2 = mean(child_age, na.rm = TRUE)
  ) %>% 
  ungroup() 

Methods

Participants

#distinct IDs at each timepoint 
cd_clean %>% 
  group_by(timepoint) %>% 
  summarise(
    n_distinct(ELS_ID)
  )

# distinct IDs at either timepoint  
cd_clean %>% 
  summarise(n_distinct(ELS_ID))

distinct_IDs_all <- 
  cd_clean %>% 
  distinct(ELS_ID) %>% 
  mutate(
    included = 1
  )

all_els_demo <-
  read_csv(els_file) %>% 
  select(ELS_ID, sumsev_threat_t1, sumsev_threat_t2, sumsev_type_t1) %>% 
  left_join(
    read_csv(demo_file, na = c("999", "888")) %>% 
      filter(timepoint == "T1"),
    by = "ELS_ID"
  ) %>% 
  left_join(distinct_IDs_all, by = "ELS_ID") %>% 
  mutate(
    included = if_else(
      is.na(included), 0, included
    )
  ) %>% 
  mutate_at(
    vars(bmi, inr_t1),
    as.double
  ) %>% 
  left_join(
    read_csv(tanner_file) %>% 
      filter(timepoint == "T1", early_puberty == TRUE) %>% 
      distinct(ELS_ID, tanner_av, .keep_all = TRUE),
    by = "ELS_ID"
  ) %>% 
  distinct(ELS_ID, sumsev_threat_t1, .keep_all = TRUE)

Parsed with column specification:
cols(
  .default = col_double(),
  class = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
  ELS_ID = col_double(),
  session_type = col_character(),
  timepoint = col_character(),
  bmi = col_character(),
  race_pr = col_character(),
  eth_pr = col_character(),
  parent_ed_t1 = col_character(),
  inr_t1 = col_character(),
  birth_date = col_character()
)
NAs introduced by coercionNAs introduced by coercionParsed with column specification:
cols(
  ELS_ID = col_double(),
  timepoint = col_character(),
  tanner_date = col_date(format = ""),
  tanner_female_pubhair = col_double(),
  tanner_female_breasts = col_double(),
  tanner_male_pubhair = col_double(),
  tanner_male_testes = col_double(),
  early_puberty = col_logical(),
  tanner_pb = col_double(),
  tanner_bt = col_double(),
  tanner_av = col_double()
)

all_els_demo %>% 
  group_by(included) %>% 
  summarise_at(
    vars(
      sumsev_threat_t1,
      sumsev_threat_t2,
      inr_t1,
      bmi,
      tanner_av
    ),
    funs(mean(., na.rm = TRUE))
  )

  
t.test(sumsev_threat_t1 ~ included, data = all_els_demo)

    Welch Two Sample t-test

data:  sumsev_threat_t1 by included
t = 0.8457, df = 40.906, p-value = 0.4026
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.6589956  1.6084338
sample estimates:
mean in group 0 mean in group 1 
       2.500000        2.025281 

t.test(sumsev_threat_t2 ~ included, data = all_els_demo)

    Welch Two Sample t-test

data:  sumsev_threat_t2 by included
t = -0.52268, df = 19.918, p-value = 0.607
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.8939999  0.5358219
sample estimates:
mean in group 0 mean in group 1 
      0.8571429       1.0362319 

t.test(bmi ~ included, data = all_els_demo)

    Welch Two Sample t-test

data:  bmi by included
t = 0.44888, df = 24.843, p-value = 0.6574
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.643051  2.558480
sample estimates:
mean in group 0 mean in group 1 
       19.13036        18.67264 

t.test(inr_t1 ~ included, data = all_els_demo)

    Welch Two Sample t-test

data:  inr_t1 by included
t = -1.6818, df = 37.965, p-value = 0.1008
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.46582485  0.04306378
sample estimates:
mean in group 0 mean in group 1 
       1.095276        1.306657 

t.test(tanner_av ~ included, data = all_els_demo)

    Welch Two Sample t-test

data:  tanner_av by included
t = -0.57163, df = 27.736, p-value = 0.5722
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.3483218  0.1963810
sample estimates:
mean in group 0 mean in group 1 
       1.863636        1.939607 

cd_wf_clean %>% 
  ggplot(aes(interval_yr)) +
  geom_histogram() +
  scale_x_continuous(breaks = seq.int(1, 3, .25)) +
  labs(
    x = "T1 to T2 interval (years)"
  )
ggsave("~/Box/lucy_king_files/ELS/cort_dhea/interval_dist.png", width = 7, height = 5)

cd_clean %>% 
  ggplot(aes(tanner_av)) +
  geom_histogram(binwidth = .5) +
  facet_grid(.~timepoint)

cd_clean %>% 
  ggplot(aes(tanner_av, fill = timepoint)) +
  geom_density(alpha = 1/2) +
  labs(
    fill = "",
    x = "Tanner stage"
  )
ggsave("~/Box/lucy_king_files/ELS/cort_dhea/tanner_dist.png", width = 7, height = 5)

cd_wf_clean %>% 
  ggplot(aes(tan_diff)) +
  geom_histogram(binwidth = .5) +
  scale_y_continuous(breaks = seq.int(0, 35, 2)) +
  labs(
    x = "Change in Tanner stage from T1 to T2"
  )
ggsave("~/Box/lucy_king_files/ELS/cort_dhea/tanner_change.png", width = 7, height = 5)

cd_clean %>% 
  filter(tan_diff == 0) %>% 
  ggplot(aes(timepoint, dhea_clean, color = factor(ELS_ID))) +
  geom_point() +
  geom_path(aes(group = ELS_ID)) 

  
cd_clean %>% 
  filter(tan_diff == 0) %>% 
  ggplot(aes(timepoint, test_clean, color = factor(ELS_ID))) +
  geom_point() +
  geom_path(aes(group = ELS_ID))

cd_wf_clean %>% 
  summarise(
    mean(tan_diff, na.rm = TRUE),
    sd(tan_diff, na.rm = TRUE)
  )

cd_descriptives <- 
  cd_clean %>% 
  filter(timepoint == "T1") %>% 
  group_by(male) %>% 
  summarise_at(
    vars(
      child_age, 
      tanner_av,
      bmi, 
      cort_win,
      dhea_win,
      test_win,
      time_collection,
      inr_t1,
      sumsev_threat_t1,
      sumsev_threat_t2
    ),
    funs(sum(!is.na(.)))
  ) %>% 
  gather(variable, n, child_age:sumsev_threat_t2) %>% 
  rename(n_T1 = n) %>% 
  left_join(
    cd_clean %>% 
      filter(timepoint == "T2") %>% 
      group_by(male) %>% 
      summarise_at(
        vars(
          child_age, 
          tanner_av,
          bmi, 
          cort_win,
          dhea_win,
          test_win,
          cort_raw,
          time_collection,
          inr_t1,
          sumsev_threat_t1,
          sumsev_threat_t2
        ),
        funs(sum(!is.na(.)))
      ) %>% 
      gather(variable, n, child_age:sumsev_threat_t2) %>% 
      rename(n_T2 = n),
    by = c("variable", "male")
  ) %>% 
  left_join(
    cd_clean %>% 
      filter(timepoint == "T1") %>% 
      group_by(male) %>% 
      summarise_at(
        vars(
          child_age, 
          tanner_av,
          bmi, 
          cort_win,
          dhea_win,
          test_win,
          time_collection,
          inr_t1,
          sumsev_threat_t1,
          sumsev_threat_t2
        ),
        funs(mean),
        na.rm = TRUE
      ) %>% 
      gather(variable, Mean, child_age:sumsev_threat_t2) %>% 
      rename(Mean_T1 = Mean),
    by = c("variable", "male")
  ) %>% 
  left_join(
    cd_clean %>% 
      filter(timepoint == "T2") %>% 
      group_by(male) %>% 
      summarise_at(
        vars(
          child_age, 
          tanner_av,
          bmi, 
          cort_win,
          dhea_win,
          test_win,
          time_collection,
          inr_t1,
          sumsev_threat_t1,
          sumsev_threat_t2
        ),
        funs(mean),
        na.rm = TRUE
      ) %>% 
      gather(variable, Mean, child_age:sumsev_threat_t2) %>% 
      rename(Mean_T2 = Mean),
    by = c("variable", "male")
  ) %>% 
  left_join(
    cd_clean %>% 
      group_by(male) %>% 
      filter(timepoint == "T1") %>% 
      summarise_at(
        vars(
          child_age, 
          tanner_av,
          bmi, 
          cort_win,
          dhea_win,
          test_win,
          time_collection,
          inr_t1,
          sumsev_threat_t1,
          sumsev_threat_t2
        ),
        funs(sd),
        na.rm = TRUE
      ) %>% 
      gather(variable, SD, child_age:sumsev_threat_t2) %>% 
    rename(SD_T1 = SD),
    by = c("variable", "male") 
  ) %>% 
  left_join(
    cd_clean %>% 
      group_by(male) %>% 
      filter(timepoint == "T2") %>% 
      summarise_at(
        vars(
          child_age, 
          tanner_av,
          bmi, 
          cort_win,
          dhea_win,
          test_win,
          time_collection,
          inr_t1,
          sumsev_threat_t1,
          sumsev_threat_t2
        ),
        funs(sd),
        na.rm = TRUE
      ) %>% 
      gather(variable, SD, child_age:sumsev_threat_t2) %>% 
      rename(SD_T2 = SD),
    by = c("variable", "male")
  ) %>% 
  mutate_at(
    vars(starts_with("Mean"), starts_with("SD")),
    funs(round), digits = 2
  ) %>% 
  unite("T1 Mean (SD)", Mean_T1, SD_T1) %>% 
  unite("T2 Mean (SD)", Mean_T2, SD_T2)  %>% 
  mutate_at(
    vars(`T1 Mean (SD)`, `T2 Mean (SD)`),
    funs(
      str_replace(
        .,
        "_", " ("
      )
    )
  ) %>% 
  mutate_at(
    vars(`T1 Mean (SD)`, `T2 Mean (SD)`),
    funs(
      paste0(. , ")")
    )
  ) %>% 
  select(
    Measure = variable,
    Sex = male,
    everything()
  ) %>% 
  mutate(
    Sex = recode(
      Sex,
      "0" = "Female",
      "1" = "Male"
    )
  )

race_ethn_descriptives <- 
  cd_clean %>% 
  filter(timepoint == "T1") %>% 
  group_by(male) %>% 
  count(
    race_pr
  ) %>% 
  mutate(
    race_pr = case_when(
      race_pr == 1 ~ "Caucasian",
      race_pr == 2 ~ "African American",
      race_pr == 3 ~ "Hispanic",
      race_pr == 4 ~ "Asian",
      race_pr == 5 ~ "Biracial",
      race_pr == 6 ~ "Other",
      TRUE ~ "Not reported"
    )
  ) %>% 
  rename(measure = race_pr) %>% 
  arrange(desc(n)) %>% 
  ungroup() %>% 
  rename(
    Measure = measure,
    n_T1 = n,
    Sex = male
  ) %>% 
  mutate(
    Sex = recode(
      Sex,
      "0" = "Female",
      "1" = "Male"
    )
  ) %>% 
  select(
    Measure,
    everything()
  )
cd_descriptives <-
  cd_descriptives %>% 
  bind_rows(race_ethn_descriptives)

cd_descriptives %>% 
  knitr::kable()

Measure	Sex	n_T1	n_T2	T1 Mean (SD)	T2 Mean (SD)
child_age	Female	94	76	11.07 (0.95)	12.96 (0.96)
child_age	Male	77	71	11.87 (0.97)	13.82 (0.98)
tanner_av	Female	94	75	2 (0.68)	3.42 (0.79)
tanner_av	Male	77	63	1.85 (0.59)	3.48 (0.92)
bmi	Female	92	62	18.31 (4.05)	20.15 (4.86)
bmi	Male	74	56	19.04 (3.38)	20.78 (3.94)
cort_win	Female	81	72	0.28 (0.18)	0.26 (0.16)
cort_win	Male	70	68	0.2 (0.11)	0.29 (0.13)
dhea_win	Female	94	74	155.81 (128.12)	168.31 (130.36)
dhea_win	Male	76	70	187.24 (330.35)	181.31 (118.9)
test_win	Female	93	76	53.03 (20.75)	61.6 (28.73)
test_win	Male	75	71	61.9 (33.77)	134.41 (76.75)
time_collection	Female	85	71	7.3 (1.2)	8.05 (1.67)
time_collection	Male	67	65	7.64 (1.29)	8.18 (1.71)
inr_t1	Female	84	68	1.32 (0.53)	1.28 (0.54)
inr_t1	Male	71	63	1.33 (0.52)	1.36 (0.51)
sumsev_threat_t1	Female	94	76	2.03 (2.13)	1.98 (2.08)
sumsev_threat_t1	Male	77	71	1.97 (1.98)	1.92 (1.87)
sumsev_threat_t2	Female	73	72	1.13 (1.99)	1.1 (1.96)
sumsev_threat_t2	Male	60	60	0.88 (1.44)	0.88 (1.43)
Caucasian	Female	46	NA	NA	NA
Caucasian	Male	35	NA	NA	NA
Biracial	Male	19	NA	NA	NA
Biracial	Female	16	NA	NA	NA
Asian	Female	15	NA	NA	NA
Hispanic	Female	8	NA	NA	NA
Asian	Male	7	NA	NA	NA
Other	Male	6	NA	NA	NA
African American	Male	5	NA	NA	NA
African American	Female	4	NA	NA	NA
Other	Female	4	NA	NA	NA
Hispanic	Male	3	NA	NA	NA
Not reported	Male	2	NA	NA	NA
Not reported	Female	1	NA	NA	NA

cd_clean %>% 
  filter(timepoint == "T1") %>% 
  count(male) %>% 
  mutate(prop = round(n / sum(n), 2))

cd_clean %>% 
  distinct(ELS_ID, inr_t1) %>% 
  count(inr_t1 < 1) %>% 
  mutate(
    prop = n / sum(n)
  )

Procedure

interval_summary <- 
  cd_wf_clean %>% 
  summarise(
    mean_interval = mean(interval_yr, na.rm = TRUE),
    sd_interval = sd(interval_yr, na.rm = TRUE),
    min_interval = min(interval_yr, na.rm = TRUE),
    max_interval = max(interval_yr, na.rm = TRUE)
  )
interval_summary

Measures

cd_clean %>% 
  filter(timepoint == "T2") %>% 
  summarise(
    min_tanner = min(tanner_av, na.rm = TRUE),
    max_tanner = max(tanner_av, na.rm = TRUE)
  )

cd_clean %>% 
  filter(timepoint == "T1") %>% 
  count(is.na(tanner_av))

cd_clean %>% 
  filter(timepoint == "T2") %>% 
  count(is.na(tanner_av))

included_IDs <- 
  cd_clean %>% 
  distinct(ELS_ID) %>% 
  pull()

read_csv(els_file) %>% 
  filter(ELS_ID %in% included_IDs) %>% 
  select(
    `Community verbal conflict ` = comm_arg_t1,
    `Community instability` = comm_inst_t1,
    `Community violence` = comm_phys_t1,
    Bullying = bully_t1,
    `Threats of domestic violence` = thrt_domvio_t1,
    `Threats of physical abuse` = thrt_phys_ab_t1,
    `Kidnapping` = kidnap_t1,
    `Family verbal conflict` =  parent_arg_t1,
    `Domestic violence` = dom_viol_t1,
    `War or terrorism` = war_terr_t1,
    `Emotional abuse` = emo_ab_t1,
    `Physical abuse` = phys_ab_t1,
    `Sexual abuse` = sex_ab_t1,
    `Witness sexual abuse` = wit_sexab_t1, 
    `Mugging or robbery` = mug_rob_t1
  ) %>% 
  gather(type, endorsed) %>% 
  count(type, endorsed) %>% 
  group_by(type) %>% 
  mutate(
    `Percent endorsed` = round(((n / sum(n)) * 100), 0)
  ) %>% 
  rename(Type = type) %>% 
  filter(endorsed == 1) %>% 
  arrange(desc(`Percent endorsed`)) %>% 
  select(-endorsed) %>% 
  knitr::kable()

Parsed with column specification:
cols(
  .default = col_double(),
  class = col_character()
)
See spec(...) for full column specifications.

Type	n	Percent endorsed
Family verbal conflict	73	41
Bullying	60	34
Community violence	21	12
Community instability	17	10
Domestic violence	14	8
Community verbal conflict	9	5
Emotional abuse	9	5
Physical abuse	8	4
Mugging or robbery	5	3
War or terrorism	5	3
Sexual abuse	3	2
Threats of domestic violence	3	2
Kidnapping	2	1
Threats of physical abuse	2	1
Witness sexual abuse	1	1

read_csv(severity_file) %>% 
  filter(ELS_ID %in% included_IDs) %>% 
  filter(timepoint == "T1") %>% 
  group_by(event, timepoint) %>% 
  summarise_at(
    vars(rating),
    funs(mean, sd)
  ) %>% 
  arrange(timepoint)

Parsed with column specification:
cols(
  ELS_ID = col_double(),
  timepoint = col_character(),
  event = col_character(),
  rating = col_double()
)

read_csv(severity_file) %>% 
  filter(ELS_ID %in% included_IDs) %>% 
  filter(timepoint == "T2") %>% 
  group_by(event, timepoint) %>% 
  summarise_at(
    vars(rating),
    funs(sum(!is.na(.)), (sum(!is.na(.))/143)*100,  mean, sd)
  ) %>% 
  arrange(timepoint)

Parsed with column specification:
cols(
  ELS_ID = col_double(),
  timepoint = col_character(),
  event = col_character(),
  rating = col_double()
)

Results

cd_clean %>% 
  group_by(timepoint) %>% 
  count(medication_binary) %>%
  mutate(prop = round(n / sum(n), 2))

Factor `medication_binary` contains implicit NA, consider using `forcats::fct_explicit_na`

cd_clean %>% 
  group_by(timepoint) %>% 
  count(corticosteroid) %>%
  mutate(prop = round(n / sum(n), 2))

cd_clean %>% 
  group_by(timepoint) %>% 
  count(birth_control) %>%
  mutate(prop = round(n / sum(n), 2))

cd_wf_clean %>% 
  select(
    sumsev_threat_t1,
    sumsev_threat_t2,
    inr_t1,
    tanner_t1,
    tanner_t2,
    bmi_t1,
    bmi_t2,
    child_age_t1,
    interval_yr,
    cort_clean_t1,
    cort_clean_t2,
    dhea_clean_t1,
    dhea_clean_t2,
    test_clean_t1,
    test_clean_t2
  ) %>% 
  correlate() %>% 
  fashion()

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

cd_wf_clean %>% 
  select(
    tanner_t1,
    tanner_t2,
    cort_clean_t1,
    cort_clean_t2,
    dhea_clean_t1,
    dhea_clean_t2,
    test_clean_t1,
    test_clean_t2
  ) %>% 
  correlate() %>% 
  fashion()

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

# DHEA ~ tanner
cor.test(cd_wf_clean$dhea_clean_t1, cd_wf_clean$tanner_t1)

    Pearson's product-moment correlation

data:  cd_wf_clean$dhea_clean_t1 and cd_wf_clean$tanner_t1
t = 1.4678, df = 168, p-value = 0.144
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.03864267  0.25865946
sample estimates:
      cor 
0.1125261 

cor.test(cd_wf_clean$dhea_clean_t2, cd_wf_clean$tanner_t2)

    Pearson's product-moment correlation

data:  cd_wf_clean$dhea_clean_t2 and cd_wf_clean$tanner_t2
t = 4.6752, df = 134, p-value = 0.000007074
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2200213 0.5106283
sample estimates:
      cor 
0.3744845 

# testosterone ~ tanner
cor.test(cd_wf_clean$test_clean_t1, cd_wf_clean$tanner_t1)

    Pearson's product-moment correlation

data:  cd_wf_clean$test_clean_t1 and cd_wf_clean$tanner_t1
t = 2.4491, df = 166, p-value = 0.01536
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.03636359 0.32885655
sample estimates:
      cor 
0.1867452 

cor.test(cd_wf_clean$test_clean_t2, cd_wf_clean$tanner_t2)

    Pearson's product-moment correlation

data:  cd_wf_clean$test_clean_t2 and cd_wf_clean$tanner_t2
t = 4.9688, df = 136, p-value = 0.000001989
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2406354 0.5247136
sample estimates:
      cor 
0.3919782 

# cortisol ~ tanner
cor.test(cd_wf_clean$cort_clean_t1, cd_wf_clean$tanner_t1)

    Pearson's product-moment correlation

data:  cd_wf_clean$cort_clean_t1 and cd_wf_clean$tanner_t1
t = 0.12029, df = 149, p-value = 0.9044
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1501105  0.1693162
sample estimates:
        cor 
0.009854267 

cor.test(cd_wf_clean$cort_clean_t2, cd_wf_clean$tanner_t2)

    Pearson's product-moment correlation

data:  cd_wf_clean$cort_clean_t2 and cd_wf_clean$tanner_t2
t = 2.3013, df = 130, p-value = 0.02296
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.02792312 0.35666605
sample estimates:
      cor 
0.1978516 

Sample characteristics

t.test(cd_wf_clean$child_age_t1, cd_wf_clean$child_age_t2, paired = TRUE)

    Paired t-test

data:  cd_wf_clean$child_age_t1 and cd_wf_clean$child_age_t2
t = -65.273, df = 141, p-value < 0.00000000000000022
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.044965 -1.924735
sample estimates:
mean of the differences 
               -1.98485 

dhea_mod_tan <-
  lmer(
    scale(dhea_clean) ~
      scale(tan_per_cent) +
      scale(tanner_av_t1t2) +
      (1|ELS_ID),
    data = cd_clean
  )
summary(dhea_mod_tan)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(dhea_clean) ~ scale(tan_per_cent) + scale(tanner_av_t1t2) +      (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 692

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5052 -0.4724  0.0067  0.4427  3.7280 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.4190   0.6473  
 Residual             0.4661   0.6827  
Number of obs: 264, groups:  ELS_ID, 133

Fixed effects:
                       Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)            -0.04149    0.07106 129.77280  -0.584 0.560278    
scale(tan_per_cent)     0.15783    0.04240 130.67196   3.723 0.000292 ***
scale(tanner_av_t1t2)   0.26917    0.07667 129.17396   3.511 0.000616 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__)
scl(tn_pr_) -0.003       
scl(tn__12) -0.159 -0.003

# sex differences: stress
t_test_stress_T1 <- t.test(cd_clean$sumsev_threat_t1 ~ cd_clean$male)
t_test_stress_T2 <- t.test(cd_clean$sumsev_threat_t2 ~ cd_clean$male)

#sex differences interval
t_test_interval <- t.test(cd_clean$interval_yr ~ cd_clean$male)

# sex differences: tanner
tanner_sex_mod <-
  lmer(
    tanner_av ~
      male * timepoint +
      (1|ELS_ID),
    data = cd_clean
  )
summary(tanner_sex_mod)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: tanner_av ~ male * timepoint + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 666.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6432 -0.4808  0.1306  0.5102  1.7782 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.2547   0.5046  
 Residual             0.2914   0.5399  
Number of obs: 309, groups:  ELS_ID, 176

Fixed effects:
                   Estimate Std. Error        df t value            Pr(>|t|)    
(Intercept)         1.98640    0.07597 264.83370  26.149 <0.0000000000000002 ***
male1              -0.12659    0.11324 264.41760  -1.118              0.2647    
timepointT2         1.41525    0.08651 155.68645  16.360 <0.0000000000000002 ***
male1:timepointT2   0.21826    0.12811 154.34916   1.704              0.0904 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) male1  tmpnT2
male1       -0.671              
timepointT2 -0.482  0.323       
ml1:tmpntT2  0.325 -0.483 -0.675

# sex differences: BMI
bmi_sex_mod <-
  lmer(
    bmi ~
      male * timepoint +
      (1|ELS_ID),
    data = cd_clean
  )
summary(bmi_sex_mod)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: bmi ~ male * timepoint + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 1492.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.77290 -0.34487 -0.07338  0.33392  2.92849 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 12.538   3.541   
 Residual              3.662   1.914   
Number of obs: 284, groups:  ELS_ID, 174

Fixed effects:
                  Estimate Std. Error       df t value             Pr(>|t|)    
(Intercept)        18.3217     0.4144 205.7578  44.212 < 0.0000000000000002 ***
male1               0.7932     0.6197 206.4266   1.280                0.202    
timepointT2         1.8735     0.3456 122.2739   5.421          0.000000303 ***
male1:timepointT2  -0.2640     0.5038 121.2448  -0.524                0.601    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) male1  tmpnT2
male1       -0.669              
timepointT2 -0.296  0.198       
ml1:tmpntT2  0.203 -0.305 -0.686

dhea_sex_mod <-
  lmer(
    dhea_clean ~
      male * timepoint +
      (1|ELS_ID),
    data = cd_clean
  )
summary(dhea_sex_mod)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: dhea_clean ~ male * timepoint + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 779.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9416 -0.4327 -0.0045  0.4631  3.7827 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.4257   0.6525  
 Residual             0.3685   0.6070  
Number of obs: 314, groups:  ELS_ID, 176

Fixed effects:
                  Estimate Std. Error        df t value             Pr(>|t|)    
(Intercept)        4.61537    0.10150 252.58631  45.474 < 0.0000000000000002 ***
male1              0.08284    0.13675 250.37416   0.606             0.545209    
timepoint1         0.38343    0.10323 140.45332   3.714             0.000293 ***
male1:timepoint1  -0.20481    0.14232 142.56740  -1.439             0.152311    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) male1  tmpnt1
male1       -0.742              
timepoint1  -0.478  0.355       
mal1:tmpnt1  0.347 -0.462 -0.725

# sex differences: cortisol
contrasts(cd_clean$timepoint) = c(0, 1) #T1 is baseline
contrasts(cd_clean$male) = c(0, 1) #Female is baseline
## timepoint____________________________________________________________________
cort_sex_mod_1 <-
  lmer(
    cort_clean ~
      male * timepoint +
      (1|ELS_ID),
    data = cd_clean
  )
summary(cort_sex_mod_1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ male * timepoint + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 547.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9426 -0.5423  0.0838  0.5875  1.8270 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.09253  0.3042  
 Residual             0.28803  0.5367  
Number of obs: 291, groups:  ELS_ID, 172

Fixed effects:
                  Estimate Std. Error        df t value             Pr(>|t|)    
(Intercept)       -1.47752    0.06829 279.89743 -21.635 < 0.0000000000000002 ***
male1             -0.28269    0.10032 279.59495  -2.818              0.00518 ** 
timepoint1        -0.02370    0.08929 152.49862  -0.265              0.79103    
male1:timepoint1   0.42361    0.12895 146.24041   3.285              0.00128 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) male1  tmpnt1
male1       -0.681              
timepoint1  -0.608  0.414       
mal1:tmpnt1  0.421 -0.615 -0.692

std_beta(cort_sex_mod_1)

contrasts(cd_clean$male) = c(1, 0) #Male is baseline
cort_sex_mod_2 <-
  lmer(
    cort_clean ~
      male * timepoint +
      (1|ELS_ID),
    data = cd_clean
  )
summary(cort_sex_mod_2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ male * timepoint + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 547.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9426 -0.5423  0.0838  0.5875  1.8270 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.09253  0.3042  
 Residual             0.28803  0.5367  
Number of obs: 291, groups:  ELS_ID, 172

Fixed effects:
                  Estimate Std. Error        df t value             Pr(>|t|)    
(Intercept)       -1.76020    0.07349 279.32906 -23.952 < 0.0000000000000002 ***
male1              0.28269    0.10032 279.59495   2.818              0.00518 ** 
timepoint1         0.39991    0.09303 140.69846   4.299            0.0000319 ***
male1:timepoint1  -0.42361    0.12895 146.24041  -3.285              0.00128 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) male1  tmpnt1
male1       -0.733              
timepoint1  -0.622  0.456       
mal1:tmpnt1  0.449 -0.615 -0.721

std_beta(cort_sex_mod_2)

# tanner________________________________________________________________________
contrasts(cd_clean$male) = c(0, 1)
cort_mod_tan_sex_girls <-
  lmer(
    scale(cort_clean) ~
      scale(tan_per_cent) * male +
      scale(tanner_av_t1t2) * male +
      (1|ELS_ID),
    data = cd_clean
  )
summary(cort_mod_tan_sex_girls)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(tan_per_cent) * male + scale(tanner_av_t1t2) *      male + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 689.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5084 -0.5703  0.1231  0.5715  1.8270 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.2247   0.474   
 Residual             0.7431   0.862   
Number of obs: 244, groups:  ELS_ID, 132

Fixed effects:
                              Estimate Std. Error         df t value Pr(>|t|)  
(Intercept)                   0.088905   0.097201 127.353804   0.915    0.362  
scale(tan_per_cent)           0.003315   0.083387 127.286290   0.040    0.968  
male1                        -0.194846   0.141134 125.077363  -1.381    0.170  
scale(tanner_av_t1t2)        -0.007914   0.104962 126.533452  -0.075    0.940  
scale(tan_per_cent):male1     0.279496   0.112070 122.994663   2.494    0.014 *
male1:scale(tanner_av_t1t2)   0.137882   0.153320 125.918785   0.899    0.370  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) male1  s(__12 s(__):
scl(tn_pr_) -0.060                            
male1       -0.689  0.041                     
scl(tn__12) -0.235  0.019  0.162              
scl(tn__):1  0.045 -0.744 -0.051 -0.014       
ml1:s(__12)  0.161 -0.013 -0.147 -0.685  0.011

#confint.merMod(cort_mod_tan_sex_girls, method = "boot")
contrasts(cd_clean$male) = c(1, 0)
cort_mod_tan_sex_boys <-
  lmer(
    scale(cort_clean) ~
      scale(tan_per_cent) * male +
      scale(tanner_av_t1t2) * male +
      (1|ELS_ID),
    data = cd_clean
  )
summary(cort_mod_tan_sex_boys)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(tan_per_cent) * male + scale(tanner_av_t1t2) *      male + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 689.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.5084 -0.5703  0.1231  0.5715  1.8270 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.2247   0.474   
 Residual             0.7431   0.862   
Number of obs: 244, groups:  ELS_ID, 132

Fixed effects:
                             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                  -0.10594    0.10233 123.05634  -1.035  0.30255    
scale(tan_per_cent)           0.28281    0.07488 117.86612   3.777  0.00025 ***
male1                         0.19485    0.14113 125.07736   1.381  0.16987    
scale(tanner_av_t1t2)         0.12997    0.11176 125.37894   1.163  0.24707    
scale(tan_per_cent):male1    -0.27950    0.11207 122.99466  -2.494  0.01396 *  
male1:scale(tanner_av_t1t2)  -0.13788    0.15332 125.91878  -0.899  0.37021    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) male1  s(__12 s(__):
scl(tn_pr_) -0.042                            
male1       -0.725  0.030                     
scl(tn__12) -0.069  0.004  0.050              
scl(tn__):1  0.028 -0.668 -0.051 -0.002       
ml1:s(__12)  0.050 -0.003 -0.147 -0.729  0.011

confint.merMod(cort_mod_tan_sex_boys, method = "boot")

Computing bootstrap confidence intervals ...

summary(test_mod_tan_girls)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(test_clean) ~ scale(tan_per_cent) * male + scale(tanner_av_t1t2) *      male + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 608.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1741 -0.5236  0.0051  0.5947  2.6971 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1982   0.4452  
 Residual             0.3901   0.6246  
Number of obs: 264, groups:  ELS_ID, 133

Fixed effects:
                             Estimate Std. Error        df t value     Pr(>|t|)    
(Intercept)                  -0.32166    0.07599 127.98902  -4.233 0.0000436645 ***
scale(tan_per_cent)           0.13739    0.05612 129.31908   2.448      0.01570 *  
male1                         0.66686    0.11103 128.95858   6.006 0.0000000181 ***
scale(tanner_av_t1t2)         0.13913    0.08234 127.98902   1.690      0.09351 .  
scale(tan_per_cent):male1     0.43047    0.07759 130.87971   5.548 0.0000001539 ***
male1:scale(tanner_av_t1t2)   0.34960    0.11974 128.43805   2.920      0.00414 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) male1  s(__12 s(__):
scl(tn_pr_)  0.000                            
male1       -0.684  0.000                     
scl(tn__12) -0.233  0.000  0.159              
scl(tn__):1  0.000 -0.723 -0.010  0.000       
ml1:s(__12)  0.160  0.000 -0.147 -0.688 -0.002

sex_diff_tp_plot <- 
  cd_clean %>% 
  select(timepoint, male, cort_clean:test_clean) %>% 
  gather(hormone, value, cort_clean:test_clean) %>% 
  mutate(
    hormone = recode(
      hormone,
      cort_clean = "log Cortisol\n(µg/dL)",
      dhea_clean = "log DHEA\n(pg/mL)",
      test_clean = "log Testosterone\n(pg/mL)"
    )
  ) %>% 
  ggplot(aes(timepoint, value, fill = male)) +
  geom_violin(draw_quantiles = .5, alpha = 3/4) +
  scale_y_continuous(breaks = seq.int(-4, 8, 1), expand = c(0, 1)) +
  scale_fill_manual(
    values = c("darkred", "darkblue"),
    labels = c("Girls", "Boys")
  ) +
  theme_apa(box = TRUE) +
  labs(
    x = "Panel A: Time-point of assessment",
    y = NULL,
    fill = NULL
  ) +
  facet_wrap(.~hormone, scales = "free") 
sex_diff_tp_plot
ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/cort_dhea_sync/sexdiff_hormones.png",
  width = 7,
  height = 5
  )

sex_diff_tan_plot <-
  cd_clean %>% 
  select(tan_per_cent, male, cort_clean:test_clean) %>% 
  gather(hormone, value, cort_clean:test_clean) %>% 
  mutate(
    hormone = recode(
      hormone,
      cort_clean = "log Cortisol\n(µg/dL)",
      dhea_clean = "log DHEA\n(pg/mL)",
      test_clean = "log Testosterone\n(pg/mL)"
    )
  ) %>% 
  ggplot(aes(tan_per_cent, value, color = male)) +
  geom_point(alpha = 1/2, size = 1.5) +
  geom_smooth(method = "lm", se = FALSE, size = 1.5) +
  scale_color_manual(
    values = c("darkred", "darkblue"),
    labels = c("Girls", "Boys")
  ) +
  theme_apa(box = TRUE) +
  labs(
    color = NULL,
    x = expression("Panel B: "*Delta*" Pubertal stage (person-mean-centered)"),
    y = NULL
  ) +
  facet_wrap(.~hormone, scales = "free") 
sex_diff_tan_plot
ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/sexdiff_hormones_tan.png",
  width = 7,
  height = 5
  )

plot_grid(sex_diff_tp_plot, sex_diff_tan_plot, rows = 2)

Argument 'rows' is deprecated. Use 'nrow' instead.Removed 34 rows containing non-finite values (stat_ydensity).Removed 182 rows containing non-finite values (stat_smooth).Removed 182 rows containing missing values (geom_point).

ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/sexdiff_hormones_grid.png",
  width = 7,
  height = 7
  )

Primary analyses

Associations between DHEA and cortisol in earlier versus later puberty: Longitudinal change and between-person Tanner stage

summary(mod1_base)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) +  
    scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 643.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.48820 -0.54027  0.07244  0.52780  2.16167 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1522   0.3901  
 Residual             0.6340   0.7962  
Number of obs: 243, groups:  ELS_ID, 132

Fixed effects:
                                         Estimate Std. Error        df t value         Pr(>|t|)    
(Intercept)                              -0.03150    0.06401 130.55905  -0.492           0.6235    
scale(dhea_clean)                         0.51806    0.06545 225.75770   7.915 0.00000000000011 ***
scale(tan_per_cent)                       0.06189    0.05458 133.96739   1.134           0.2588    
scale(tanner_av_t1t2)                    -0.05293    0.06943 128.96335  -0.762           0.4473    
scale(dhea_clean):scale(tan_per_cent)     0.14579    0.06168 178.72244   2.364           0.0192 *  
scale(dhea_clean):scale(tanner_av_t1t2)   0.03832    0.06862 236.88212   0.558           0.5771    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(_) sc(__) s(__12 s(_):(__)
scl(dh_cln) -0.038                               
scl(tn_pr_) -0.004 -0.214                        
scl(tn__12) -0.136 -0.206  0.065                 
scl(_):(__) -0.118  0.193 -0.248 -0.109          
s(_):(__12) -0.190  0.124 -0.134 -0.020  0.085   

# collection time_________________________________
cd_clean_time <-
  cd_clean %>% 
  filter(!is.na(timecoll_av_t1t2))
mod1_time <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)
mod1_base_time <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(timecoll_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)
anova(mod1_time, mod1_base_time)

refitting model(s) with ML (instead of REML)

Data: cd_clean_time
Models:
mod1_time: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_time:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_time: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_time:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(timecoll_av_t1t2) + 
mod1_base_time:     (1 | ELS_ID)
               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_time       8 630.70 658.58 -307.35   614.70                         
mod1_base_time  9 632.55 663.91 -307.27   614.55 0.1481      1     0.7004

# interaction with change in tanner
mod1_base_timeX1 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * scale(timecoll_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)
anova(mod1_time, mod1_base_timeX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean_time
Models:
mod1_time: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_time:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_timeX1: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_timeX1:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_base_timeX1:     scale(timecoll_av_t1t2) + (1 | ELS_ID)
                 Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_time         8 630.70 658.58 -307.35   614.70                         
mod1_base_timeX1 10 633.65 668.50 -306.83   613.65 1.0445      2     0.5932

# interaction with mean tanner
mod1_base_timeX2 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * scale(timecoll_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)
anova(mod1_time, mod1_base_timeX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean_time
Models:
mod1_time: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_time:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_timeX2: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_timeX2:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_base_timeX2:     scale(timecoll_av_t1t2) + (1 | ELS_ID)
                 Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_time         8 630.70 658.58 -307.35   614.70                         
mod1_base_timeX2 10 634.47 669.32 -307.24   614.47 0.2249      2     0.8936

# T1 to T2 interval________________________________
cd_clean_interval <-
  cd_clean %>% 
  filter(!is.na(interval_yr))
mod1_int <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_interval
)
mod1_base_int <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(interval_yr) +
    (1 | ELS_ID),
  data = cd_clean_interval
)
anova(mod1_int, mod1_base_int)

refitting model(s) with ML (instead of REML)

Data: cd_clean_interval
Models:
mod1_int: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_int:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_int: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_int:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(interval_yr) + 
mod1_base_int:     (1 | ELS_ID)
              Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_int       8 631.16 658.97 -307.58   615.16                         
mod1_base_int  9 633.03 664.32 -307.52   615.03 0.1231      1     0.7256

# interaction with change in tanner
mod1_base_intX1 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * scale(interval_yr) +
    (1 | ELS_ID),
  data = cd_clean_interval
)
anova(mod1_int, mod1_base_intX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean_interval
Models:
mod1_int: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_int:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_intX1: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_intX1:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_base_intX1:     scale(interval_yr) + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_int         8 631.16 658.97 -307.58   615.16                         
mod1_base_intX1 10 631.24 666.01 -305.62   611.24 3.9133      2     0.1413

# interaction with mean tanner
mod1_base_intX2 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2)  * scale(interval_yr) +
    (1 | ELS_ID),
  data = cd_clean_interval
)
anova(mod1_int, mod1_base_intX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean_interval
Models:
mod1_int: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_int:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_intX2: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_intX2:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_base_intX2:     scale(interval_yr) + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_int         8 631.16 658.97 -307.58   615.16                         
mod1_base_intX2 10 635.02 669.78 -307.51   615.02 0.1403      2     0.9323

# medication________________________________
mod1_med <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    medication_binary +
    (1 | ELS_ID),
  data = cd_clean
)
    
anova(mod1_base, mod1_med)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_med: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_med:     scale(dhea_clean) * scale(tanner_av_t1t2) + medication_binary + 
mod1_med:     (1 | ELS_ID)
          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_base  8 636.83 664.77 -310.41   620.83                         
mod1_med   9 638.83 670.27 -310.41   620.83 0.0005      1     0.9818

# interaction with change in tanner
mod1_medX1 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * medication_binary +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_base, mod1_medX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_medX1: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_medX1:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_medX1:     medication_binary + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_base   8 636.83 664.77 -310.41   620.83                         
mod1_medX1 10 639.06 673.99 -309.53   619.06 1.7674      2     0.4133

# interaction with mean tanner
mod1_medX2 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * medication_binary +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_base, mod1_medX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_medX2: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_medX2:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_medX2:     medication_binary + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_base   8 636.83 664.77 -310.41   620.83                         
mod1_medX2 10 640.82 675.75 -310.41   620.82 0.0085      2     0.9957

# sex________________________________________
mod1_sex <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    male +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_base, mod1_sex)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_sex: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_sex:     scale(dhea_clean) * scale(tanner_av_t1t2) + male + (1 | ELS_ID)
          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_base  8 636.83 664.77 -310.41   620.83                         
mod1_sex   9 636.68 668.12 -309.34   618.68 2.1497      1     0.1426

# interaction with change in tanner
mod1_sexX1 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * male +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_base, mod1_sexX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_sexX1: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_sexX1:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_sexX1:     male + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
mod1_base   8 636.83 664.77 -310.41   620.83                           
mod1_sexX1 10 633.15 668.08 -306.57   613.15 7.6787      2    0.02151 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# interaction with mean tanner
mod1_sexX2 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * male +
    scale(tanner_av_t1t2) * male +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_sexX1, mod1_sexX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_sexX1: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_sexX1:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_sexX1:     male + (1 | ELS_ID)
mod1_sexX2: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_sexX2:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_sexX2:     male + scale(tanner_av_t1t2) * male + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_sexX1 10 633.15 668.08 -306.57   613.15                         
mod1_sexX2 11 634.93 673.35 -306.46   612.93 0.2226      1     0.6371

# check results hold when controlling for interaction with sex -- YES 
summary(mod1_sexX1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) +  
    scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) *      male + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 643.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.59708 -0.51026  0.05349  0.58514  2.02338 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1524   0.3903  
 Residual             0.6165   0.7852  
Number of obs: 243, groups:  ELS_ID, 132

Fixed effects:
                                         Estimate Std. Error        df t value          Pr(>|t|)    
(Intercept)                              -0.03936    0.06364 129.30059  -0.619            0.5373    
scale(dhea_clean)                         0.50977    0.06485 224.08648   7.861 0.000000000000159 ***
scale(tan_per_cent)                       0.05082    0.05400 132.68998   0.941            0.3484    
scale(tanner_av_t1t2)                    -0.06188    0.06906 127.70075  -0.896            0.3719    
male1                                    -0.09787    0.06193 125.31705  -1.580            0.1166    
scale(dhea_clean):scale(tan_per_cent)     0.14980    0.06113 175.43166   2.450            0.0153 *  
scale(dhea_clean):scale(tanner_av_t1t2)   0.04882    0.06801 234.93905   0.718            0.4735    
scale(tan_per_cent):male1                 0.12017    0.05143 122.36902   2.337            0.0211 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(_) sc(__) s(__12 male1  s(_):(__) s(_):(__1
scl(dh_cln) -0.038                                                
scl(tn_pr_) -0.003 -0.210                                         
scl(tn__12) -0.127 -0.205  0.066                                  
male1        0.087 -0.003  0.026  0.091                           
scl(_):(__) -0.124  0.194 -0.246 -0.116 -0.084                    
s(_):(__12) -0.194  0.122 -0.137 -0.025 -0.069  0.089             
scl(tn__):1 -0.001 -0.059 -0.075  0.001 -0.059 -0.023     0.027   

# bmi T1______________________________________
cd_clean_bmi <-
  cd_clean %>% 
  filter(!is.na(bmi_t1))
mod1_bmi <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
mod1_base_bmi <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(bmi_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
    
anova(mod1_bmi, mod1_base_bmi)

refitting model(s) with ML (instead of REML)

Data: cd_clean_bmi
Models:
mod1_bmi: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_bmi:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_bmi: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_bmi:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(bmi_av_t1t2) + 
mod1_base_bmi:     (1 | ELS_ID)
              Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_bmi       8 615.29 643.03 -299.64   599.29                         
mod1_base_bmi  9 615.21 646.43 -298.61   597.21 2.0753      1     0.1497

# interaction with change in tanner
mod1_base_bmiX1 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * scale(bmi_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
    
anova(mod1_bmi, mod1_base_bmiX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean_bmi
Models:
mod1_bmi: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_bmi:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_bmiX1: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_bmiX1:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_base_bmiX1:     scale(bmi_av_t1t2) + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_bmi         8 615.29 643.03 -299.64   599.29                         
mod1_base_bmiX1 10 615.58 650.26 -297.79   595.58 3.7121      2     0.1563

# interaction with mean tanner
mod1_base_bmiX2 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * scale(bmi_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
    
anova(mod1_bmi, mod1_base_bmiX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean_bmi
Models:
mod1_bmi: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_bmi:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_base_bmiX2: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base_bmiX2:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_base_bmiX2:     scale(bmi_av_t1t2) + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_bmi         8 615.29 643.03 -299.64   599.29                         
mod1_base_bmiX2 10 617.17 651.85 -298.58   597.17 2.1203      2     0.3464

# age________________________________________
mod1_age <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(age_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_base, mod1_age)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_age: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_age:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(age_av_t1t2) + 
mod1_age:     (1 | ELS_ID)
          Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_base  8 636.83 664.77 -310.41   620.83                         
mod1_age   9 638.58 670.02 -310.29   620.58 0.2463      1     0.6197

# interaction with change in tanner
mod1_ageX1 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * scale(age_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_base, mod1_ageX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_ageX1: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_ageX1:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_ageX1:     scale(age_av_t1t2) + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_base   8 636.83 664.77 -310.41   620.83                         
mod1_ageX1 10 638.90 673.83 -309.45   618.90 1.9252      2     0.3819

# interaction with mean tanner
mod1_ageX2 <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
   scale(tanner_av_t1t2) * scale(age_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_base, mod1_ageX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_base: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_base:     scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_ageX2: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) + 
mod1_ageX2:     scale(dhea_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_ageX2:     scale(age_av_t1t2) + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_base   8 636.83 664.77 -310.41   620.83                         
mod1_ageX2 10 640.41 675.34 -310.20   620.41 0.4229      2     0.8094

cd_clean_nosteroid <-
  cd_clean %>% 
  filter(corticosteroid == 0)
mod1_nosteroid <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_nosteroid
)
summary(mod1_nosteroid)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_clean) * scale(tan_per_cent) +      scale(dhea_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
   Data: cd_clean_nosteroid

REML criterion at convergence: 606.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.51331 -0.54110  0.07302  0.49285  2.19475 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1508   0.3883  
 Residual             0.6433   0.8021  
Number of obs: 228, groups:  ELS_ID, 129

Fixed effects:
                                         Estimate Std. Error        df t value         Pr(>|t|)    
(Intercept)                              -0.02427    0.06608 125.24502  -0.367           0.7140    
scale(dhea_clean)                         0.49534    0.06751 209.02477   7.337 0.00000000000475 ***
scale(tan_per_cent)                       0.08991    0.05751 131.93692   1.563           0.1204    
scale(tanner_av_t1t2)                    -0.05783    0.07160 122.85670  -0.808           0.4208    
scale(dhea_clean):scale(tan_per_cent)     0.13037    0.06348 167.32342   2.054           0.0416 *  
scale(dhea_clean):scale(tanner_av_t1t2)   0.01792    0.07107 221.68331   0.252           0.8012    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(_) sc(__) s(__12 s(_):(__)
scl(dh_cln) -0.029                               
scl(tn_pr_) -0.006 -0.224                        
scl(tn__12) -0.134 -0.207  0.069                 
scl(_):(__) -0.125  0.191 -0.262 -0.110          
s(_):(__12) -0.190  0.119 -0.134 -0.024  0.094   

#confint.merMod(t1t2_mod1_nosteroid, method = "boot")

confint.merMod(mod1_fin_lo, method = "boot")

Computing bootstrap confidence intervals ...

7 message(s): boundary (singular) fit: see ?isSingular

                                                                                2.5 %     97.5 %
.sig01                                                                     0.15977381 0.57881527
.sigma                                                                     0.67629121 0.89065871
(Intercept)                                                               -0.24993800 0.08177118
scale(dhea_clean)                                                          0.23527107 0.54620644
I(scale(tan_per_cent) + sd(tan_per_cent, na.rm = TRUE))                   -0.04266982 0.16417884
scale(tanner_av_t1t2)                                                     -0.19070069 0.08938766
scale(dhea_clean):I(scale(tan_per_cent) + sd(tan_per_cent, na.rm = TRUE))  0.01538456 0.25599204
scale(dhea_clean):scale(tanner_av_t1t2)                                   -0.09930373 0.17211155

cd_clean_boys <-
  cd_clean %>% 
  filter(male == 1)
mod1_boys_hi <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * I(scale(tan_per_cent) - sd(tan_per_cent, na.rm = TRUE)) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_boys
)
summary(mod1_boys_hi)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_clean) * I(scale(tan_per_cent) -  
    sd(tan_per_cent, na.rm = TRUE)) + scale(dhea_clean) * scale(tanner_av_t1t2) +      (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 281.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2239 -0.5668  0.0776  0.5247  1.9008 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1855   0.4307  
 Residual             0.4832   0.6951  
Number of obs: 112, groups:  ELS_ID, 60

Fixed effects:
                                                                           Estimate Std. Error        df t value
(Intercept)                                                                 0.07820    0.11499  98.99346   0.680
scale(dhea_clean)                                                           0.83992    0.14651 104.26599   5.733
I(scale(tan_per_cent) - sd(tan_per_cent, na.rm = TRUE))                     0.15474    0.07274  62.48155   2.127
scale(tanner_av_t1t2)                                                      -0.05758    0.09941  59.48193  -0.579
scale(dhea_clean):I(scale(tan_per_cent) - sd(tan_per_cent, na.rm = TRUE))   0.24777    0.09204  80.82245   2.692
scale(dhea_clean):scale(tanner_av_t1t2)                                     0.02603    0.09427 105.94630   0.276
                                                                              Pr(>|t|)    
(Intercept)                                                                    0.49804    
scale(dhea_clean)                                                         0.0000000968 ***
I(scale(tan_per_cent) - sd(tan_per_cent, na.rm = TRUE))                        0.03735 *  
scale(tanner_av_t1t2)                                                          0.56465    
scale(dhea_clean):I(scale(tan_per_cent) - sd(tan_per_cent, na.rm = TRUE))      0.00863 ** 
scale(dhea_clean):scale(tanner_av_t1t2)                                        0.78298    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(_) I-sn=T s(__12 s-sn=T
scl(dh_cln) -0.366                            
I((__)-sn=T  0.609 -0.404                     
scl(tn__12)  0.000 -0.266  0.104              
s(_):I-sn=T -0.346  0.794 -0.320 -0.148       
s(_):(__12) -0.314  0.328 -0.204 -0.097  0.207

#confint.merMod(mod1_boys_hi, method = "boot")
mod1_boys_hi_int <- summary(mod1_boys_hi)$coefficients[1]
mod1_boys_hi_slp <- summary(mod1_boys_hi)$coefficients[2]
mod1_boys_lo <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * I(scale(tan_per_cent) + sd(tan_per_cent, na.rm = TRUE)) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_boys
)
summary(mod1_boys_lo)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_clean) * I(scale(tan_per_cent) +  
    sd(tan_per_cent, na.rm = TRUE)) + scale(dhea_clean) * scale(tanner_av_t1t2) +      (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 281.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.2239 -0.5668  0.0776  0.5247  1.9008 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1855   0.4307  
 Residual             0.4832   0.6951  
Number of obs: 112, groups:  ELS_ID, 60

Fixed effects:
                                                                           Estimate Std. Error        df t value Pr(>|t|)
(Intercept)                                                                -0.20612    0.11119  98.51126  -1.854 0.066770
scale(dhea_clean)                                                           0.38466    0.10347 100.34534   3.718 0.000331
I(scale(tan_per_cent) + sd(tan_per_cent, na.rm = TRUE))                     0.15474    0.07274  62.48155   2.127 0.037352
scale(tanner_av_t1t2)                                                      -0.05758    0.09941  59.48193  -0.579 0.564649
scale(dhea_clean):I(scale(tan_per_cent) + sd(tan_per_cent, na.rm = TRUE))   0.24777    0.09204  80.82245   2.692 0.008628
scale(dhea_clean):scale(tanner_av_t1t2)                                     0.02603    0.09427 105.94630   0.276 0.782976
                                                                             
(Intercept)                                                               .  
scale(dhea_clean)                                                         ***
I(scale(tan_per_cent) + sd(tan_per_cent, na.rm = TRUE))                   *  
scale(tanner_av_t1t2)                                                        
scale(dhea_clean):I(scale(tan_per_cent) + sd(tan_per_cent, na.rm = TRUE)) ** 
scale(dhea_clean):scale(tanner_av_t1t2)                                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(_) I+sn=T s(__12 s+sn=T
scl(dh_cln)  0.106                            
I((__)+sn=T -0.572 -0.048                     
scl(tn__12) -0.125 -0.135  0.104              
s(_):I+sn=T  0.028 -0.510 -0.320 -0.148       
s(_):(__12) -0.079  0.126 -0.204 -0.097  0.207

#confint.merMod(mod1_boys_lo, method = "boot")
mod1_boys_lo_int <- summary(mod1_boys_lo)$coefficients[1]
mod1_boys_lo_slp <- summary(mod1_boys_lo)$coefficients[2]

cd_clean_girls <-
  cd_clean %>% 
  filter(male == 0)

mod1_girls <- lmer(
  scale(cort_clean) ~ 
    scale(dhea_clean) * scale(tan_per_cent) +
    scale(dhea_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_girls
)

summary(mod1_girls)
confint.merMod(mod1_girls, method = "boot")

Associations between testosterone and cortisol in earlier versus later puberty: Longitudinal change and between-person Tanner stage

confint.merMod(mod1_test_base, method = "boot")

Computing bootstrap confidence intervals ...

2 warning(s): Model failed to converge with max|grad| = 0.002233 (tol = 0.002, component 1) (and others)

                                             2.5 %      97.5 %
.sig01                                   0.2907855  0.63856874
.sigma                                   0.6650252  0.86483130
(Intercept)                             -0.1169910  0.16362691
scale(test_clean)                        0.3308260  0.60461693
scale(tan_per_cent)                     -0.1095842  0.10470578
scale(tanner_av_t1t2)                   -0.2031724  0.11849515
scale(test_clean):scale(tan_per_cent)   -0.1020184  0.13631871
scale(test_clean):scale(tanner_av_t1t2) -0.3359531 -0.02538013

# collection time_________________________________
cd_clean_time <-
  cd_clean %>% 
  filter(!is.na(timecoll_av_t1t2))

mod1_test_time <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)

mod1_test_base_time <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(timecoll_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)

anova(mod1_test_time, mod1_test_base_time)

# interaction with change in tanner

mod1_test_base_timeX1 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    tan_per_cent * scale(timecoll_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)
anova(mod1_test_time, mod1_test_base_timeX1)

# interaction with mean tanner

mod1_test_base_timeX2 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * scale(timecoll_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_time
)
anova(mod1_test_time, mod1_test_base_timeX2)

# T1 to T2 interval________________________________
cd_clean_interval <-
  cd_clean %>% 
  filter(!is.na(interval_yr))

mod1_test_int <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_interval
)

mod1_test_base_int <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(interval_yr) +
    (1 | ELS_ID),
  data = cd_clean_interval
)

anova(mod1_test_int, mod1_test_base_int)

# interaction with change in tanner

mod1_test_base_intX1 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * interval_yr +
    (1 | ELS_ID),
  data = cd_clean_interval
)
anova(mod1_test_int, mod1_test_base_intX1)

# interaction with mean tanner
mod1_test_base_intX2 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * interval_yr +
    (1 | ELS_ID),
  data = cd_clean_interval
)
anova(mod1_test_int, mod1_test_base_intX2)

# medication________________________________
mod1_test_med <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    medication_binary +
    (1 | ELS_ID),
  data = cd_clean
)
    
anova(mod1_test_base, mod1_test_med)

# interaction with change in tanner
mod1_test_medX1 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * medication_binary +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_test_base, mod1_test_medX1)

# interaction with mean tanner
mod1_test_medX2 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * medication_binary +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_test_base, mod1_test_medX2)

# sex________________________________________
mod1_test_sex <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    male +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_test_base, mod1_test_sex)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_test_base: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_base:     scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_test_sex: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_sex:     scale(test_clean) * scale(tanner_av_t1t2) + male + (1 | ELS_ID)
               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
mod1_test_base  8 649.98 677.96 -316.99   633.98                             
mod1_test_sex   9 636.52 668.00 -309.26   618.52 15.461      1 0.00008424 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# interaction with change in tanner
mod1_test_sexX1 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    tan_per_cent * male +
    (1 | ELS_ID),
  data = cd_clean
)

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

anova(mod1_test_sex, mod1_test_sexX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_test_sex: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_sex:     scale(test_clean) * scale(tanner_av_t1t2) + male + (1 | ELS_ID)
mod1_test_sexX1: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_sexX1:     scale(test_clean) * scale(tanner_av_t1t2) + tan_per_cent * 
mod1_test_sexX1:     male + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_sex    9 636.52 668.00 -309.26   618.52                         
mod1_test_sexX1 10 638.44 673.41 -309.22   618.44 0.0827      1     0.7737

# interaction with mean tanner
mod1_test_sexX2 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * male +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_test_sex, mod1_test_sexX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_test_sex: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_sex:     scale(test_clean) * scale(tanner_av_t1t2) + male + (1 | ELS_ID)
mod1_test_sexX2: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_sexX2:     scale(test_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_test_sexX2:     male + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_sex    9 636.52 668.00 -309.26   618.52                         
mod1_test_sexX2 10 638.50 673.47 -309.25   618.50 0.0208      1     0.8853

#check results hold when including sex -- NO, become marginal
summary(mod1_test_sex)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) +  
    scale(test_clean) * scale(tanner_av_t1t2) + male + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 644.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.81073 -0.56210  0.02634  0.59866  2.24580 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1731   0.416   
 Residual             0.6053   0.778   
Number of obs: 244, groups:  ELS_ID, 132

Fixed effects:
                                          Estimate Std. Error         df t value          Pr(>|t|)    
(Intercept)                              -0.007539   0.067682 147.865583  -0.111          0.911460    
scale(test_clean)                         0.564070   0.073830 236.943831   7.640 0.000000000000535 ***
scale(tan_per_cent)                      -0.049361   0.057964 147.204972  -0.852          0.395830    
scale(tanner_av_t1t2)                    -0.108172   0.071347 132.653306  -1.516          0.131862    
male1                                    -0.275788   0.069208 147.627033  -3.985          0.000106 ***
scale(test_clean):scale(tan_per_cent)     0.044475   0.060074 194.307385   0.740          0.459985    
scale(test_clean):scale(tanner_av_t1t2)  -0.138323   0.072592 232.890814  -1.905          0.057949 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(_) sc(__) s(__12 male1  s(_):(__)
scl(tst_cl)  0.090                                      
scl(tn_pr_)  0.035 -0.389                               
scl(tn__12) -0.133 -0.266  0.139                        
male1        0.134 -0.313  0.217  0.177                 
scl(_):(__) -0.291 -0.256 -0.159 -0.032 -0.135          
s(_):(__12) -0.187 -0.079 -0.082  0.048 -0.142 -0.135   

confint.merMod(mod1_test_sex, method = "boot")

Computing bootstrap confidence intervals ...

1 message(s): boundary (singular) fit: see ?isSingular
1 warning(s): Model failed to converge with max|grad| = 0.00368574 (tol = 0.002, component 1)

                                              2.5 %       97.5 %
.sig01                                   0.21233575  0.589142862
.sigma                                   0.66161916  0.866726504
(Intercept)                             -0.14516127  0.127757273
scale(test_clean)                        0.41882938  0.704280142
scale(tan_per_cent)                     -0.16024359  0.061137833
scale(tanner_av_t1t2)                   -0.24652991  0.038841787
male1                                   -0.40961003 -0.129357160
scale(test_clean):scale(tan_per_cent)   -0.07580428  0.167571431
scale(test_clean):scale(tanner_av_t1t2) -0.28304987  0.002327858

# bmi T1______________________________________
cd_clean_bmi <-
  cd_clean %>% 
  filter(!is.na(bmi_t1))
mod1_test_bmi <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
mod1_test_base_bmi <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(bmi_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
    
anova(mod1_test_bmi, mod1_test_base_bmi)

refitting model(s) with ML (instead of REML)

Data: cd_clean_bmi
Models:
mod1_test_bmi: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_bmi:     scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_test_base_bmi: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_base_bmi:     scale(test_clean) * scale(tanner_av_t1t2) + scale(bmi_av_t1t2) + 
mod1_test_base_bmi:     (1 | ELS_ID)
                   Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_bmi       8 633.10 660.88 -308.55   617.10                         
mod1_test_base_bmi  9 634.37 665.62 -308.19   616.37 0.7289      1     0.3932

# interaction with change in tanner
mod1_test_base_bmiX1 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * scale(bmi_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
    
anova(mod1_test_bmi, mod1_test_base_bmiX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean_bmi
Models:
mod1_test_bmi: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_bmi:     scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_test_base_bmiX1: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_base_bmiX1:     scale(test_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_test_base_bmiX1:     scale(bmi_av_t1t2) + (1 | ELS_ID)
                     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_bmi         8 633.10 660.88 -308.55   617.10                         
mod1_test_base_bmiX1 10 635.28 670.00 -307.64   615.28 1.8216      2     0.4022

# interaction with mean tanner
mod1_test_base_bmiX2 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * scale(bmi_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_bmi
)
    
anova(mod1_test_bmi, mod1_test_base_bmiX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean_bmi
Models:
mod1_test_bmi: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_bmi:     scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_test_base_bmiX2: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_base_bmiX2:     scale(test_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_test_base_bmiX2:     scale(bmi_av_t1t2) + (1 | ELS_ID)
                     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_bmi         8 633.10 660.88 -308.55   617.10                         
mod1_test_base_bmiX2 10 634.11 668.84 -307.06   614.11 2.9889      2     0.2244

# age________________________________________
mod1_test_age <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(age_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_test_base, mod1_test_age)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_test_base: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_base:     scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_test_age: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_age:     scale(test_clean) * scale(tanner_av_t1t2) + scale(age_av_t1t2) + 
mod1_test_age:     (1 | ELS_ID)
               Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_base  8 649.98 677.96 -316.99   633.98                         
mod1_test_age   9 651.11 682.59 -316.56   633.11 0.8726      1     0.3502

# interaction with change in tanner
mod1_test_ageX1 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tan_per_cent) * scale(age_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_test_base, mod1_test_ageX1)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_test_base: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_base:     scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_test_ageX1: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_ageX1:     scale(test_clean) * scale(tanner_av_t1t2) + scale(tan_per_cent) * 
mod1_test_ageX1:     scale(age_av_t1t2) + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_base   8 649.98 677.96 -316.99   633.98                         
mod1_test_ageX1 10 651.14 686.11 -315.57   631.14 2.8426      2     0.2414

# interaction with mean tanner
mod1_test_ageX2 <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    scale(tanner_av_t1t2) * scale(age_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean
)
anova(mod1_test_base, mod1_test_ageX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod1_test_base: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_base:     scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
mod1_test_ageX2: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) + 
mod1_test_ageX2:     scale(test_clean) * scale(tanner_av_t1t2) + scale(tanner_av_t1t2) * 
mod1_test_ageX2:     scale(age_av_t1t2) + (1 | ELS_ID)
                Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod1_test_base   8 649.98 677.96 -316.99   633.98                         
mod1_test_ageX2 10 653.11 688.08 -316.56   633.11 0.8734      2     0.6462

cd_clean_nosteroid <-
  cd_clean %>% 
  filter(corticosteroid == 0)
contrasts(cd_clean_nosteroid$timepoint) = c(0, 1) # T1 is baseline
mod1_test_nosteroid <- lmer(
  scale(cort_clean) ~ 
    scale(test_clean) * scale(tan_per_cent) +
    scale(test_clean) * scale(tanner_av_t1t2) +
    (1 | ELS_ID),
  data = cd_clean_nosteroid
)
summary(mod1_test_nosteroid)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_clean) * scale(tan_per_cent) +      scale(test_clean) * scale(tanner_av_t1t2) + (1 | ELS_ID)
   Data: cd_clean_nosteroid

REML criterion at convergence: 613.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.61180 -0.55355 -0.00011  0.58209  2.28610 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.2321   0.4817  
 Residual             0.5953   0.7716  
Number of obs: 229, groups:  ELS_ID, 129

Fixed effects:
                                         Estimate Std. Error        df t value       Pr(>|t|)    
(Intercept)                               0.02825    0.07257 142.19283   0.389         0.6977    
scale(test_clean)                         0.48281    0.07497 222.91070   6.440 0.000000000723 ***
scale(tan_per_cent)                       0.01300    0.05977 144.42957   0.217         0.8282    
scale(tanner_av_t1t2)                    -0.07315    0.07585 126.38920  -0.964         0.3367    
scale(test_clean):scale(tan_per_cent)     0.02477    0.06182 178.52518   0.401         0.6891    
scale(test_clean):scale(tanner_av_t1t2)  -0.17731    0.07588 221.82489  -2.337         0.0204 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(_) sc(__) s(__12 s(_):(__)
scl(tst_cl)  0.148                               
scl(tn_pr_)  0.007 -0.370                        
scl(tn__12) -0.155 -0.224  0.116                 
scl(_):(__) -0.284 -0.308 -0.144 -0.026          
s(_):(__12) -0.169 -0.145 -0.053  0.074 -0.156   

confint.merMod(mod1_test_fin_lo, method = "boot")

Computing bootstrap confidence intervals ...

                                                                                   2.5 %      97.5 %
.sig01                                                                         0.3010394  0.65844977
.sigma                                                                         0.6714998  0.87579900
(Intercept)                                                                   -0.1185596  0.26024616
scale(test_clean)                                                              0.4092759  0.77080184
scale(tan_per_cent)                                                           -0.1125931  0.10436265
I(scale(tanner_av_t1t2) + sd(tanner_av_t1t2, na.rm = TRUE))                   -0.2151385  0.08814345
scale(test_clean):scale(tan_per_cent)                                         -0.1048126  0.13354803
scale(test_clean):I(scale(tanner_av_t1t2) + sd(tanner_av_t1t2, na.rm = TRUE)) -0.3210396 -0.03443752

confint.merMod(mod1_test_boys_lo, method = "boot")

Computing bootstrap confidence intervals ...

6 message(s): boundary (singular) fit: see ?isSingular

                                                                                   2.5 %      97.5 %
.sig01                                                                         0.1812263  0.64259391
.sigma                                                                         0.5539406  0.81500898
(Intercept)                                                                   -0.1338990  0.34770765
scale(test_clean)                                                              0.5612616  1.09073534
scale(tan_per_cent)                                                           -0.1764934  0.19186567
I(scale(tanner_av_t1t2) + sd(tanner_av_t1t2, na.rm = TRUE))                   -0.3544284  0.04892210
scale(test_clean):scale(tan_per_cent)                                         -0.1089150  0.26009367
scale(test_clean):I(scale(tanner_av_t1t2) + sd(tanner_av_t1t2, na.rm = TRUE)) -0.4672985 -0.06227382

confint.merMod(mod1_test_girls, method = "boot")

Computing bootstrap confidence intervals ...

36 message(s): boundary (singular) fit: see ?isSingular

                                             2.5 %     97.5 %
.sig01                                   0.0000000 0.65765525
.sigma                                   0.6736694 0.99478709
(Intercept)                             -0.1539868 0.19480789
scale(test_clean)                        0.1826737 0.57292116
scale(tan_per_cent)                     -0.2242877 0.08144782
scale(tanner_av_t1t2)                   -0.2631407 0.15904266
scale(test_clean):scale(tan_per_cent)   -0.1797972 0.16223477
scale(test_clean):scale(tanner_av_t1t2) -0.1925920 0.28689632

plot_mod1_test_sex <- 
  cd_clean %>% 
  filter(!is.na(tanner_av)) %>% 
  ggplot(
    aes(
      test_clean,
      cort_clean,
      color = factor(ceiling(tanner_av) >= 3)
    )
  ) + 
  geom_point(size = 1.5, alpha = 1/2) +
  geom_smooth(method = "lm", se = FALSE, size = 1.5) +
  scale_x_continuous(breaks = seq.int(1, 8, 1), expand = c(0, 1)) +
  scale_y_continuous(breaks = seq.int(-4, 1, 1), expand = c(0, 1)) +
  scale_color_manual(
    values = c("darkblue", "darkred"),
    labels = c("Pubertal stages 1-2", "Pubertal stages 3-5")
  ) +
  theme_apa(box = TRUE) +
  theme(
    legend.position = "bottom",
    legend.text = element_text(size = 12)
  ) +
  labs(
    color = NULL,
    x = "Panel B: log Testosterone (pg/mL)",
    y = "log Cortisol (µg/dL)"
  ) +
  facet_wrap(
    .~ male, 
    scales = "free",
    labeller = labeller(male = c("0" = "Girls", "1" = "Boys"))
  )
plot_mod1_test_sex
ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/cort_dhea_sync/cort_test_tan_sex.png",
  height = 4, width = 5
)

plot_grid(plot_mod1_sex_tp, plot_mod1_test_sex, rows = 2)
ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/cort_dhea_sync/cort_hormone_tan_sex_grid.png",
  height = 7, width = 5
)

Association of longitudinal change and mean DHEA with cortisol

confint.merMod(mod2_base, method = "boot")

Computing bootstrap confidence intervals ...

5 message(s): boundary (singular) fit: see ?isSingular

                           2.5 %    97.5 %
.sig01                0.15605528 0.5687590
.sigma                0.67904497 0.8818697
(Intercept)          -0.09871565 0.1141914
scale(dhea_per_cent)  0.21270267 0.4021419
scale(dhea_av_t1t2)   0.24071837 0.5223217

anova(mod2_time_fit, mod2_time)

refitting model(s) with ML (instead of REML)

Data: cd_clean_time
Models:
mod2_time_fit: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) + 
mod2_time_fit:     (1 | ELS_ID)
mod2_time: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) + 
mod2_time:     scale(timecoll_av_t1t2) + (1 | ELS_ID)
              Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod2_time_fit  5 658.23 675.90 -324.11   648.23                         
mod2_time      6 660.05 681.25 -324.03   648.05 0.1749      1     0.6758

anova(mod2_base, mod2_medX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod2_base: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) + 
mod2_base:     (1 | ELS_ID)
mod2_medX2: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) * 
mod2_medX2:     medication_binary + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod2_base   5 664.48 682.19 -327.24   654.48                         
mod2_medX2  7 667.87 692.66 -326.93   653.87 0.6143      2     0.7355

anova(mod2_base, mod2_sexX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod2_base: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) + 
mod2_base:     (1 | ELS_ID)
mod2_sexX2: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) * 
mod2_sexX2:     male + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod2_base   5 664.48 682.19 -327.24   654.48                         
mod2_sexX2  7 667.13 691.92 -326.56   653.13 1.3569      2     0.5074

anova(mod2_bmi_fit, mod2_bmiX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean_bmi
Models:
mod2_bmi_fit: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) + 
mod2_bmi_fit:     (1 | ELS_ID)
mod2_bmiX2: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) * 
mod2_bmiX2:     scale(bmi_av_t1t2) + (1 | ELS_ID)
             Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod2_bmi_fit  5 657.33 675.00 -323.67   647.33                         
mod2_bmiX2    7 659.40 684.13 -322.70   645.40 1.9339      2     0.3802

anova(mod2_base, mod2_ageX2)

refitting model(s) with ML (instead of REML)

Data: cd_clean
Models:
mod2_base: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) + 
mod2_base:     (1 | ELS_ID)
mod2_ageX2: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) * 
mod2_ageX2:     scale(age_av_t1t2) + (1 | ELS_ID)
           Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
mod2_base   5 664.48 682.19 -327.24   654.48                         
mod2_ageX2  7 667.89 692.68 -326.94   653.89 0.5963      2     0.7422

mod2_nosteroid_fit <-
 lmer(
    scale(cort_clean) ~
      scale(dhea_per_cent) +
      scale(dhea_av_t1t2) +
      (1 | ELS_ID),
    data = cd_clean_nosteroid
  )
summary(mod2_nosteroid_fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) +      (1 | ELS_ID)
   Data: cd_clean_nosteroid

REML criterion at convergence: 630.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.86180 -0.51176  0.04587  0.53945  2.11295 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1542   0.3926  
 Residual             0.6408   0.8005  
Number of obs: 240, groups:  ELS_ID, 134

Fixed effects:
                       Estimate Std. Error         df t value     Pr(>|t|)    
(Intercept)            0.006142   0.062206 129.588259   0.099        0.922    
scale(dhea_per_cent)   0.302072   0.052717 119.919753   5.730 0.0000000761 ***
scale(dhea_av_t1t2)    0.365953   0.068684 141.202168   5.328 0.0000003832 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__)
scl(dh_pr_) -0.024       
scl(dh__12) -0.008  0.018

summary(mod2_boys)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) +      (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 308.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0536 -0.4515  0.1386  0.5131  1.5122 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1729   0.4159  
 Residual             0.5272   0.7261  
Number of obs: 123, groups:  ELS_ID, 65

Fixed effects:
                      Estimate Std. Error        df t value     Pr(>|t|)    
(Intercept)           0.002491   0.083695 62.856828   0.030        0.976    
scale(dhea_per_cent)  0.407723   0.066670 59.893555   6.116 0.0000000792 ***
scale(dhea_av_t1t2)   0.382037   0.092651 71.120602   4.123 0.0000998153 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__)
scl(dh_pr_) -0.025       
scl(dh__12)  0.012 -0.025

summary(mod2_girls)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_per_cent) + scale(dhea_av_t1t2) +      (1 | ELS_ID)
   Data: cd_clean_girls

REML criterion at convergence: 360.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.53601 -0.57673  0.00757  0.61165  1.88020 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1544   0.393   
 Residual             0.7174   0.847   
Number of obs: 132, groups:  ELS_ID, 72

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           0.009414   0.087542 67.257263   0.108 0.914684    
scale(dhea_per_cent)  0.192638   0.075635 67.141084   2.547 0.013165 *  
scale(dhea_av_t1t2)   0.354865   0.098347 74.595800   3.608 0.000555 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__)
scl(dh_pr_) -0.029       
scl(dh__12) -0.028  0.067

DHEA/cort: Effects of ELS threat severity

confint.merMod(mod2_threat, method = "boot")

Computing bootstrap confidence intervals ...

                                                   2.5 %      97.5 %
.sig01                                        0.24469321  0.59114865
.sigma                                        0.66745456  0.85811543
(Intercept)                                  -0.12704630  0.11913730
scale(dhea_per_cent)                          0.22204888  0.42586189
scale(sumsev_threat_t1)                      -0.17067610  0.06726760
scale(dhea_av_t1t2)                           0.24381887  0.51815200
scale(dhea_per_cent):scale(sumsev_threat_t1) -0.24636331 -0.02804225
scale(sumsev_threat_t1):scale(dhea_av_t1t2)  -0.06632965  0.17856716

confint.merMod(mod2_threat, method = "boot")

Computing bootstrap confidence intervals ...

1 message(s): boundary (singular) fit: see ?isSingular

                                                   2.5 %      97.5 %
.sig01                                        0.23531218  0.59957630
.sigma                                        0.65986511  0.87075875
(Intercept)                                  -0.13194838  0.10838763
scale(dhea_per_cent)                          0.22455186  0.42611698
scale(sumsev_threat_t1)                      -0.17926115  0.05825371
scale(dhea_av_t1t2)                           0.23757418  0.52242607
scale(dhea_per_cent):scale(sumsev_threat_t1) -0.24824742 -0.02573026
scale(sumsev_threat_t1):scale(dhea_av_t1t2)  -0.06571214  0.19000949

summary(mod2_threat_t2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_per_cent) * scale(sumsev_threat_t2) +  
    scale(dhea_av_t1t2) * scale(sumsev_threat_t2) + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 620.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.86349 -0.49875  0.00442  0.56704  2.02437 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1569   0.3961  
 Residual             0.6587   0.8116  
Number of obs: 231, groups:  ELS_ID, 124

Fixed effects:
                                              Estimate Std. Error        df t value     Pr(>|t|)    
(Intercept)                                   -0.02267    0.06446 117.88728  -0.352        0.726    
scale(dhea_per_cent)                           0.33326    0.05400 113.57548   6.172 0.0000000107 ***
scale(sumsev_threat_t2)                       -0.04405    0.06617 122.61207  -0.666        0.507    
scale(dhea_av_t1t2)                            0.39122    0.07527 128.99647   5.197 0.0000007698 ***
scale(dhea_per_cent):scale(sumsev_threat_t2)  -0.03363    0.04752 114.24846  -0.708        0.481    
scale(sumsev_threat_t2):scale(dhea_av_t1t2)    0.08827    0.10195 127.04314   0.866        0.388    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) sc(__2) s(__12 s(__):
scl(dh_pr_) -0.023                             
scl(sms__2)  0.012  0.010                      
scl(dh__12)  0.012  0.032  0.036               
s(__):(__2)  0.016 -0.145  0.055   0.045       
s(__2):(__1  0.004  0.037  0.096   0.373  0.072

mod2_threat_boys <-
 lmer(
    scale(cort_clean) ~
      scale(dhea_per_cent) * scale(sumsev_threat_t1) +
      scale(dhea_av_t1t2) * scale(sumsev_threat_t1) +
      (1 | ELS_ID),
    data = cd_clean_boys
  )
summary(mod2_threat_boys)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_per_cent) * scale(sumsev_threat_t1) +  
    scale(dhea_av_t1t2) * scale(sumsev_threat_t1) + (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 309.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1811 -0.4430  0.1299  0.5833  1.8452 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1460   0.3822  
 Residual             0.5212   0.7219  
Number of obs: 123, groups:  ELS_ID, 65

Fixed effects:
                                             Estimate Std. Error       df t value     Pr(>|t|)    
(Intercept)                                  -0.03907    0.08268 60.58499  -0.473       0.6382    
scale(dhea_per_cent)                          0.41245    0.06631 58.84922   6.220 0.0000000561 ***
scale(sumsev_threat_t1)                      -0.15267    0.08734 60.88996  -1.748       0.0855 .  
scale(dhea_av_t1t2)                           0.44766    0.09324 66.91835   4.801 0.0000092389 ***
scale(dhea_per_cent):scale(sumsev_threat_t1) -0.07261    0.07449 59.81648  -0.975       0.3336    
scale(sumsev_threat_t1):scale(dhea_av_t1t2)   0.17004    0.08040 63.18517   2.115       0.0384 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) sc(__1) s(__12 s(__):
scl(dh_pr_) -0.026                             
scl(sms__1)  0.075 -0.025                      
scl(dh__12) -0.034 -0.018 -0.238               
s(__):(__1) -0.027 -0.023 -0.037   0.003       
s(__1):(__1 -0.197  0.002 -0.052   0.149 -0.028

#confint.merMod(mod2_threat_boys, method = "boot")
# higher threat severity in boys________________________________________________
mod2_threat_boys_hi <-
   lmer(
    cort_clean ~
      scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(dhea_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_boys
  )
summary(mod2_threat_boys_hi)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) + scale(dhea_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) -  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 205.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1811 -0.4430  0.1299  0.5833  1.8452 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.05962  0.2442  
 Residual             0.21276  0.4613  
Number of obs: 123, groups:  ELS_ID, 65

Fixed effects:
                                                                                                                   Estimate
(Intercept)                                                                                                        -1.68940
scale(dhea_per_cent, scale = FALSE)                                                                                 0.42669
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     -0.05071
scale(dhea_av_t1t2, scale = FALSE)                                                                                  0.49421
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) -0.04740
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)   0.07073
                                                                                                                   Std. Error
(Intercept)                                                                                                           0.07965
scale(dhea_per_cent, scale = FALSE)                                                                                   0.12379
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                        0.02901
scale(dhea_av_t1t2, scale = FALSE)                                                                                    0.10551
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))    0.04862
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)     0.03344
                                                                                                                         df
(Intercept)                                                                                                        62.94528
scale(dhea_per_cent, scale = FALSE)                                                                                60.60006
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     60.88996
scale(dhea_av_t1t2, scale = FALSE)                                                                                 61.60929
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) 59.81648
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  63.18517
                                                                                                                   t value
(Intercept)                                                                                                        -21.210
scale(dhea_per_cent, scale = FALSE)                                                                                  3.447
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                      -1.748
scale(dhea_av_t1t2, scale = FALSE)                                                                                   4.684
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))  -0.975
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)    2.115
                                                                                                                               Pr(>|t|)
(Intercept)                                                                                                        < 0.0000000000000002
scale(dhea_per_cent, scale = FALSE)                                                                                             0.00104
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                                  0.08550
scale(dhea_av_t1t2, scale = FALSE)                                                                                            0.0000159
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))              0.33359
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)               0.03839
                                                                                                                      
(Intercept)                                                                                                        ***
scale(dhea_per_cent, scale = FALSE)                                                                                ** 
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     .  
scale(dhea_av_t1t2, scale = FALSE)                                                                                 ***
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))    
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F-sn=T s(__1s=F ss=Fs=F-sn=T
s(__,s=FALS  -0.056                                         
I((s=F-sn=T   0.750 -0.045                                  
s(__12,s=FA  -0.236 -0.019   -0.200                         
ss=Fs=F-sn=T -0.044  0.740   -0.037    -0.015               
Is=F-sn=Ts=F -0.167 -0.020   -0.052     0.715   -0.028      

mod2_threat_boys_hi_int <- summary(mod2_threat_boys_hi)$coefficients[1]
mod2_threat_boys_hi_slp <- summary(mod2_threat_boys_hi)$coefficients[2]
# lower threat severity in boys________________________________________________
mod2_threat_boys_lo <-
   lmer(
    cort_clean ~
      scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(dhea_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_boys
  )
summary(mod2_threat_boys_lo)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) + scale(dhea_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) +  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 205.3

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1811 -0.4430  0.1299  0.5833  1.8452 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.05962  0.2442  
 Residual             0.21276  0.4613  
Number of obs: 123, groups:  ELS_ID, 65

Fixed effects:
                                                                                                                   Estimate
(Intercept)                                                                                                        -1.49431
scale(dhea_per_cent, scale = FALSE)                                                                                 0.60904
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     -0.05071
scale(dhea_av_t1t2, scale = FALSE)                                                                                  0.22212
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) -0.04740
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)   0.07073
                                                                                                                   Std. Error
(Intercept)                                                                                                           0.07392
scale(dhea_per_cent, scale = FALSE)                                                                                   0.12663
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                        0.02901
scale(dhea_av_t1t2, scale = FALSE)                                                                                    0.09097
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))    0.04862
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)     0.03344
                                                                                                                         df
(Intercept)                                                                                                        58.28468
scale(dhea_per_cent, scale = FALSE)                                                                                58.24995
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     60.88996
scale(dhea_av_t1t2, scale = FALSE)                                                                                 70.57110
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) 59.81648
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  63.18517
                                                                                                                   t value
(Intercept)                                                                                                        -20.214
scale(dhea_per_cent, scale = FALSE)                                                                                  4.810
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                      -1.748
scale(dhea_av_t1t2, scale = FALSE)                                                                                   2.442
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))  -0.975
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)    2.115
                                                                                                                               Pr(>|t|)
(Intercept)                                                                                                        < 0.0000000000000002
scale(dhea_per_cent, scale = FALSE)                                                                                           0.0000111
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                                   0.0855
scale(dhea_av_t1t2, scale = FALSE)                                                                                               0.0171
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))               0.3336
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)                0.0384
                                                                                                                      
(Intercept)                                                                                                        ***
scale(dhea_per_cent, scale = FALSE)                                                                                ***
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     .  
scale(dhea_av_t1t2, scale = FALSE)                                                                                 *  
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))    
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F+sn=T s(__1s=F ss=Fs=F+sn=T
s(__,s=FALS  -0.006                                         
I((s=F+sn=T  -0.702  0.011                                  
s(__12,s=FA   0.199 -0.027   -0.158                         
ss=Fs=F+sn=T  0.008 -0.754   -0.037     0.022               
Is=F+sn=Ts=F -0.101  0.022   -0.052    -0.585   -0.028      

mod2_threat_boys_lo_int <- summary(mod2_threat_boys_lo)$coefficients[1]
mod2_threat_boys_lo_slp <- summary(mod2_threat_boys_lo)$coefficients[2]

mod2_threat_girls <-
 lmer(
    scale(cort_clean) ~
      scale(dhea_per_cent) * scale(sumsev_threat_t1) +
      scale(dhea_av_t1t2) * scale(sumsev_threat_t1) +
      (1 | ELS_ID),
    data = cd_clean_girls
  )
summary(mod2_threat_girls)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(dhea_per_cent) * scale(sumsev_threat_t1) +  
    scale(dhea_av_t1t2) * scale(sumsev_threat_t1) + (1 | ELS_ID)
   Data: cd_clean_girls

REML criterion at convergence: 364.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.38209 -0.56675 -0.01006  0.58978  1.74404 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1908   0.4368  
 Residual             0.6736   0.8208  
Number of obs: 132, groups:  ELS_ID, 72

Fixed effects:
                                              Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                                   0.006990   0.088694 65.690218   0.079 0.937422    
scale(dhea_per_cent)                          0.230829   0.075612 67.048811   3.053 0.003249 ** 
scale(sumsev_threat_t1)                      -0.007694   0.091763 64.150751  -0.084 0.933440    
scale(dhea_av_t1t2)                           0.353129   0.099552 72.715607   3.547 0.000686 ***
scale(dhea_per_cent):scale(sumsev_threat_t1) -0.176580   0.077109 63.505476  -2.290 0.025354 *  
scale(sumsev_threat_t1):scale(dhea_av_t1t2)  -0.029030   0.096546 63.028165  -0.301 0.764645    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) sc(__1) s(__12 s(__):
scl(dh_pr_) -0.033                             
scl(sms__1)  0.019  0.006                      
scl(dh__12) -0.025  0.060 -0.053               
s(__):(__1)  0.016 -0.235  0.011   0.011       
s(__1):(__1 -0.039  0.024  0.149  -0.053  0.000

#confint.merMod(mod2_threat_girls, method = "boot")
# higher threat severity in girls________________________________________________
mod2_threat_girls_hi <-
   lmer(
    cort_clean ~
      scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(dhea_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_girls
  )
summary(mod2_threat_girls_hi)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) + scale(dhea_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) -  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_girls

REML criterion at convergence: 244.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.38209 -0.56675 -0.01006  0.58978  1.74404 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.0739   0.2718  
 Residual             0.2609   0.5108  
Number of obs: 132, groups:  ELS_ID, 72

Fixed effects:
                                                                                                                    Estimate
(Intercept)                                                                                                        -1.483702
scale(dhea_per_cent, scale = FALSE)                                                                                 0.086336
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     -0.002274
scale(dhea_av_t1t2, scale = FALSE)                                                                                  0.260522
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) -0.133444
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  -0.011081
                                                                                                                   Std. Error
(Intercept)                                                                                                          0.080178
scale(dhea_per_cent, scale = FALSE)                                                                                  0.150330
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                       0.027119
scale(dhea_av_t1t2, scale = FALSE)                                                                                   0.108501
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))   0.058273
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)    0.036852
                                                                                                                          df
(Intercept)                                                                                                        64.806814
scale(dhea_per_cent, scale = FALSE)                                                                                62.719645
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     64.150751
scale(dhea_av_t1t2, scale = FALSE)                                                                                 67.309140
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) 63.505476
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  63.028165
                                                                                                                   t value
(Intercept)                                                                                                        -18.505
scale(dhea_per_cent, scale = FALSE)                                                                                  0.574
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                      -0.084
scale(dhea_av_t1t2, scale = FALSE)                                                                                   2.401
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))  -2.290
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)   -0.301
                                                                                                                              Pr(>|t|)
(Intercept)                                                                                                        <0.0000000000000002
scale(dhea_per_cent, scale = FALSE)                                                                                             0.5678
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                                  0.9334
scale(dhea_av_t1t2, scale = FALSE)                                                                                              0.0191
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))              0.0254
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)               0.7646
                                                                                                                      
(Intercept)                                                                                                        ***
scale(dhea_per_cent, scale = FALSE)                                                                                   
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                        
scale(dhea_av_t1t2, scale = FALSE)                                                                                 *  
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) *  
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F-sn=T s(__1s=F ss=Fs=F-sn=T
s(__,s=FALS  0.001                                          
I((s=F-sn=T  0.725  0.014                                   
s(__12,s=FA  0.016  0.056    0.068                          
ss=Fs=F-sn=T 0.019  0.628    0.011     0.009                
Is=F-sn=Ts=F 0.079  0.020    0.149     0.676    0.000       

mod2_threat_girls_hi_int <- summary(mod2_threat_girls_hi)$coefficients[1]
mod2_threat_girls_hi_slp <- summary(mod2_threat_girls_hi)$coefficients[2]
# lower threat severity in girls________________________________________________
mod2_threat_girls_lo <-
   lmer(
    cort_clean ~
      scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(dhea_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_girls
  )
summary(mod2_threat_girls_lo)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(dhea_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) + scale(dhea_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) +  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_girls

REML criterion at convergence: 244.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.38209 -0.56675 -0.01006  0.58978  1.74404 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.0739   0.2718  
 Residual             0.2609   0.5108  
Number of obs: 132, groups:  ELS_ID, 72

Fixed effects:
                                                                                                                    Estimate
(Intercept)                                                                                                        -1.474125
scale(dhea_per_cent, scale = FALSE)                                                                                 0.648381
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     -0.002274
scale(dhea_av_t1t2, scale = FALSE)                                                                                  0.307192
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) -0.133444
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  -0.011081
                                                                                                                   Std. Error
(Intercept)                                                                                                          0.078672
scale(dhea_per_cent, scale = FALSE)                                                                                  0.190999
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                       0.027119
scale(dhea_av_t1t2, scale = FALSE)                                                                                   0.114371
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))   0.058273
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)    0.036852
                                                                                                                          df
(Intercept)                                                                                                        64.976305
scale(dhea_per_cent, scale = FALSE)                                                                                66.815857
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     64.150751
scale(dhea_av_t1t2, scale = FALSE)                                                                                 68.362770
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) 63.505476
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)  63.028165
                                                                                                                   t value
(Intercept)                                                                                                        -18.738
scale(dhea_per_cent, scale = FALSE)                                                                                  3.395
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                      -0.084
scale(dhea_av_t1t2, scale = FALSE)                                                                                   2.686
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))  -2.290
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)   -0.301
                                                                                                                               Pr(>|t|)
(Intercept)                                                                                                        < 0.0000000000000002
scale(dhea_per_cent, scale = FALSE)                                                                                             0.00116
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                                  0.93344
scale(dhea_av_t1t2, scale = FALSE)                                                                                              0.00907
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))              0.02535
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)               0.76464
                                                                                                                      
(Intercept)                                                                                                        ***
scale(dhea_per_cent, scale = FALSE)                                                                                ** 
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                        
scale(dhea_av_t1t2, scale = FALSE)                                                                                 ** 
scale(dhea_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) *  
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(dhea_av_t1t2, scale = FALSE)     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F+sn=T s(__1s=F ss=Fs=F+sn=T
s(__,s=FALS  -0.020                                         
I((s=F+sn=T  -0.713 -0.003                                  
s(__12,s=FA   0.106  0.011   -0.138                         
ss=Fs=F+sn=T  0.003 -0.791    0.011     0.008               
Is=F+sn=Ts=F -0.135  0.015    0.149    -0.715    0.000      

mod2_threat_girls_lo_int <- summary(mod2_threat_girls_lo)$coefficients[1]
mod2_threat_girls_lo_slp <- summary(mod2_threat_girls_lo)$coefficients[2]

plot_mod2 <-
  tibble(
    Sex = c("Girls", "Boys", "Girls", "Boys"),
    Threat = c("Higher (+1SD)", "Higher (+1SD)", "Lower (-1SD)", "Lower (-1SD)"),
    Intercept = c(
      mod2_threat_girls_hi_int, 
      mod2_threat_boys_hi_int, 
      mod2_threat_girls_lo_int, 
      mod2_threat_boys_lo_int
    ),
    Slope = c(
      mod2_threat_girls_hi_slp, 
      mod2_threat_boys_hi_slp, 
      mod2_threat_girls_lo_slp, 
      mod2_threat_boys_lo_slp
    )
  )
# produce a plot for legend only_______________________________________________
plot_mod2_threat_legend <- 
  cd_clean_boys %>% 
  ggplot(aes(dhea_per_cent, cort_clean)) +
  geom_point(alpha = 1/2) +
  geom_abline(
    data = plot_mod2 %>% filter(Sex == "Boys"),
    aes(
      intercept = Intercept,
      slope = Slope,
      color = Threat
    ),
    size = 2
  ) +
  scale_x_continuous(breaks = seq.int(-4, 4, 1)) +
  scale_color_manual(
    values = c("darkred", "darkblue")
  ) +
  expand_limits(y = c(-1, -4)) +
  expand_limits(x = c(-2, 2)) +
  theme_apa(base_size = 20, box = TRUE) +
  theme(
    plot.title = element_text(hjust = .5)
  ) +
  labs(
    title = "Boys",
    shape = NULL,
    size = NULL,
    color = "Threat severity",
    x = NULL,
    y = "log Cortisol (µg/dL)"
  )  
# plot for boys________________________________________________________________
plot_mod2_threat_boys <- 
  cd_clean_boys %>% 
  ggplot(aes(dhea_per_cent, cort_clean)) +
  geom_point(alpha = 1/2, aes(size = sumsev_threat_t1, show.legend = FALSE)) +
  geom_abline(
    data = plot_mod2 %>% filter(Sex == "Boys"),
    aes(
      intercept = Intercept,
      slope = Slope,
      color = Threat
    ),
    size = 2
  ) +
  scale_x_continuous(breaks = seq.int(-4, 4, 1)) +
  scale_color_manual(
    values = c("darkred", "darkblue")
  ) +
  expand_limits(y = c(-1, -4)) +
  expand_limits(x = c(-2, 2)) +
  theme_apa(base_size = 20, box = TRUE) +
  theme(
    plot.title = element_text(hjust = .5)
  ) +
  labs(
    title = "Boys",
    shape = NULL,
    size = NULL,
    color = "Threat severity",
    x = NULL,
    y = NULL
  )  

Ignoring unknown aesthetics: show.legend

ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/dhea_cort_threat_boys.png",
  height = 6, width = 10
)
# plot for girls________________________________________________________________
plot_mod2_threat_girls <- 
  cd_clean_girls %>% 
  ggplot(aes(dhea_per_cent, cort_clean)) +
  geom_point(alpha = 1/2, aes(size = sumsev_threat_t1), show.legend = FALSE) +
  geom_abline(
    data = plot_mod2 %>% filter(Sex == "Girls"),
    aes(
      intercept = Intercept,
      slope = Slope,
      color = Threat
    ),
    size = 2
  ) +
  scale_x_continuous(breaks = seq.int(-4, 4, 1)) +
  scale_color_manual(
    values = c("darkred", "darkblue")
  ) +
  expand_limits(y = c(-1, -4)) +
  expand_limits(x = c(-2, 2)) +
  theme_apa(base_size = 20, box = TRUE) +
  theme(
    plot.title = element_text(hjust = .5)
  ) +
  labs(
    title = "Girls",
    shape = NULL,
    size = NULL,
    color = "Threat severity",
    x = NULL,
    y = NULL
  )  
ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/dhea_cort_threat_girls.png",
  height = 6, width = 10
)
# save together________________________________________________________________
plot_mod2_threat_sex <- 
  plot_grid(
    plot_mod2_threat_girls + theme(legend.position = "none"),
    plot_mod2_threat_boys + theme(legend.position = "none"), 
    align = "hv", 
    axis = "bt", 
    rel_widths = c(1, 1)
  ) 
#create common x and y labels___________________________________________________
library(grid)
library(gridExtra)
library(gridGraphics)
y.grob_cort <- textGrob(
  "log Cortisol (µg/dL)",
  gp = gpar(col = "black", fontsize = 22), 
  rot = 90
)
x.grob_dhea <- textGrob(
  expression("Panel A: "*Delta*" log DHEA (pg/mL) (person-mean-centered)"),
  gp = gpar(col = "black", fontsize = 22)
)
#add common labels to plot______________________________________________________
plot_mod2_threat_sex <-
  grid.arrange(
  arrangeGrob(
    plot_mod2_threat_sex, 
    left = y.grob_cort, 
    bottom = x.grob_dhea
  )
)

# extract legend________________________________________________________________
threat_legend <- get_legend(
  # create some space to the left of the legend
  plot_mod2_threat_legend + 
    theme(
      legend.box.margin = margin(0, 0, 0, 12)
    )
)

Removed 25 rows containing missing values (geom_point).

# add legend back_______________________________________________________________
# the width of one plot (via rel_widths).
plot_grid(plot_mod2_threat_sex, threat_legend, rel_widths = c(1, .2))
# save plot_____________________________________________________________________
ggsave("~/Box/lucy_king_files/ELS/cort_dhea/dhea_cort_threat_sex.png", width = 12, height = 6)

Association of longitudinal change and mean testosterone with cortisol

mod2_test_base <-
 lmer(
    scale(cort_clean) ~
      scale(test_per_cent) +
      scale(test_av_t1t2) +
      (1 | ELS_ID),
    data = cd_clean
  )
summary(mod2_test_base)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_per_cent) + scale(test_av_t1t2) +      (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 676.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.70678 -0.56209  0.04319  0.56189  2.48626 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.2249   0.4742  
 Residual             0.5993   0.7742  
Number of obs: 256, groups:  ELS_ID, 137

Fixed effects:
                       Estimate Std. Error         df t value      Pr(>|t|)    
(Intercept)            0.003143   0.063433 133.358986   0.050      0.960562    
scale(test_per_cent)   0.313906   0.048834 126.044478   6.428 0.00000000242 ***
scale(test_av_t1t2)    0.278008   0.069974 136.458734   3.973      0.000114 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__)
scl(tst_p_) -0.025       
scl(ts__12) -0.005 -0.013

cd_clean_time <-
  cd_clean %>% 
  filter(!is.na(timecoll_av_t1t2))

mod2_test_time_fit <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_time
  )

mod2_test_time <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      timecoll_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_time
  )

anova(mod2_test_time_fit, mod2_test_time)

mod2_test_timeX1 <-
 lmer(
    cort_clean ~
      test_per_cent * timecoll_av_t1t2 +
      test_av_t1t2 +
      timecoll_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_time
  )

anova(mod2_test_time_fit, mod2_test_timeX1)

mod2_test_timeX2 <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 * timecoll_av_t1t2 +
      timecoll_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_time
  )

anova(mod2_test_time_fit, mod2_test_timeX2)

mod2_test_med <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      medication_binary +
      (1 | ELS_ID),
    data = cd_clean
  )

anova(mod2_test_base, mod2_test_med)

mod2_test_medX1 <-
 lmer(
    cort_clean ~
      test_per_cent * medication_binary +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean
  )

anova(mod2_test_base, mod2_test_medX1)

mod2_test_medX2 <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 * medication_binary +
      (1 | ELS_ID),
    data = cd_clean
  )

anova(mod2_test_base, mod2_test_medX2)

mod2_test_sex <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      male +
      (1 | ELS_ID),
    data = cd_clean
  )

anova(mod2_test_base, mod2_test_sex)

mod2_test_sexX1 <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean
  )
summary(mod2_test_sexX1)

anova(mod2_test_base, mod2_test_sexX1)

mod2_test_sexX2 <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 * male +
      (1 | ELS_ID),
    data = cd_clean
  )

anova(mod2_test_base, mod2_test_sexX2)

cd_clean_bmi <-
  cd_clean %>% 
  filter(!is.na(bmi_av_t1t2))

mod2_test_bmi_fit <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_bmi
  )

mod2_test_bmi <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      bmi_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_bmi
  )

anova(mod2_test_bmi_fit, mod2_test_bmi)

mod2_test_bmiX1 <-
 lmer(
    cort_clean ~
      test_per_cent * bmi_av_t1t2 +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_bmi
  )

anova(mod2_test_bmi_fit, mod2_test_bmiX1)

mod2_test_bmiX2 <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 * bmi_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_bmi
  )

anova(mod2_test_bmi_fit, mod2_test_bmiX2)

mod2_test_nosteroid_fit <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_nosteroid
  )

summary(mod2_test_nosteroid_fit)

#confint.merMod(t1t2_mod1_test_nosteroid, method = "boot")

mod2_test_boys <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_boys
  )
summary(mod2_test_boys)
#confint.merMod(mod2_test_boys, method = "boot")

mod2_test_girls <-
 lmer(
    cort_clean ~
      test_per_cent +
      test_av_t1t2 +
      (1 | ELS_ID),
    data = cd_clean_girls
  )
summary(mod2_test_girls)

confint.merMod(mod2_test_girls, method = "boot")

Testosteron/cort: Effects of ELS threat severity

summary(mod2_test_threat)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_per_cent) * scale(sumsev_threat_t1) +  
    scale(test_av_t1t2) * scale(sumsev_threat_t1) + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 681.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.78906 -0.57460  0.02752  0.59678  2.28282 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.2441   0.4941  
 Residual             0.5763   0.7592  
Number of obs: 256, groups:  ELS_ID, 137

Fixed effects:
                                               Estimate Std. Error         df t value      Pr(>|t|)    
(Intercept)                                    0.004361   0.063900 131.288257   0.068       0.94569    
scale(test_per_cent)                           0.307242   0.048004 124.817480   6.400 0.00000000284 ***
scale(sumsev_threat_t1)                       -0.023839   0.065974 129.820787  -0.361       0.71843    
scale(test_av_t1t2)                            0.284041   0.071270 134.453435   3.985       0.00011 ***
scale(test_per_cent):scale(sumsev_threat_t1)  -0.116697   0.048175 124.470920  -2.422       0.01686 *  
scale(sumsev_threat_t1):scale(test_av_t1t2)   -0.012293   0.063104 131.428195  -0.195       0.84585    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) sc(__1) s(__12 s(__):
scl(tst_p_) -0.026                             
scl(sms__1)  0.036 -0.015                      
scl(ts__12) -0.003 -0.011 -0.045               
s(__):(__1) -0.012  0.058 -0.029  -0.013       
s(__1):(__1 -0.013 -0.017  0.145  -0.153 -0.033

confint.merMod(mod2_test_threat, method = "boot")

Computing bootstrap confidence intervals ...

1 warning(s): Model failed to converge with max|grad| = 0.00253107 (tol = 0.002, component 1)

                                                  2.5 %      97.5 %
.sig01                                        0.3154895  0.66061863
.sigma                                        0.6670611  0.85778966
(Intercept)                                  -0.1261734  0.14042704
scale(test_per_cent)                          0.2164492  0.40035091
scale(sumsev_threat_t1)                      -0.1559008  0.09747200
scale(test_av_t1t2)                           0.1595182  0.43561248
scale(test_per_cent):scale(sumsev_threat_t1) -0.2123697 -0.02130992
scale(sumsev_threat_t1):scale(test_av_t1t2)  -0.1321208  0.10965877

mod2_test_threat_t2 <-
 lmer(
    scale(cort_clean) ~
      scale(test_per_cent) * scale(sumsev_threat_t2) +
      scale(test_av_t1t2) * scale(sumsev_threat_t2) +
      (1 | ELS_ID),
    data = cd_clean
  )
summary(mod2_test_threat_t2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_per_cent) * scale(sumsev_threat_t2) +  
    scale(test_av_t1t2) * scale(sumsev_threat_t2) + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 632.3

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.55994 -0.53660  0.01771  0.58563  2.46461 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.2325   0.4822  
 Residual             0.6338   0.7961  
Number of obs: 232, groups:  ELS_ID, 124

Fixed effects:
                                              Estimate Std. Error        df t value      Pr(>|t|)    
(Intercept)                                   -0.02027    0.06828 118.54793  -0.297      0.767042    
scale(test_per_cent)                           0.34323    0.05462 112.55750   6.284 0.00000000644 ***
scale(sumsev_threat_t2)                       -0.04804    0.06976 121.55874  -0.689      0.492341    
scale(test_av_t1t2)                            0.27122    0.07576 118.73172   3.580      0.000499 ***
scale(test_per_cent):scale(sumsev_threat_t2)  -0.03177    0.05442 110.63505  -0.584      0.560510    
scale(sumsev_threat_t2):scale(test_av_t1t2)    0.03813    0.09247 117.00845   0.412      0.680836    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) sc(__2) s(__12 s(__):
scl(tst_p_) -0.019                             
scl(sms__2)  0.022 -0.002                      
scl(ts__12)  0.028  0.004  0.026               
s(__):(__2) -0.002  0.018 -0.014   0.012       
s(__2):(__1  0.032  0.015  0.075  -0.005  0.011

mod2_test_threat_boys <-
 lmer(
    scale(cort_clean) ~
      scale(test_per_cent) * scale(sumsev_threat_t1) +
      scale(test_av_t1t2) * scale(sumsev_threat_t1) +
      (1 | ELS_ID),
    data = cd_clean_boys
  )
summary(mod2_test_threat_boys)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_per_cent) * scale(sumsev_threat_t1) +  
    scale(test_av_t1t2) * scale(sumsev_threat_t1) + (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 305.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.40026 -0.48992  0.08387  0.64318  1.62534 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1622   0.4027  
 Residual             0.4865   0.6975  
Number of obs: 123, groups:  ELS_ID, 64

Fixed effects:
                                             Estimate Std. Error       df t value     Pr(>|t|)    
(Intercept)                                  -0.04994    0.08277 59.67464  -0.603       0.5485    
scale(test_per_cent)                          0.41020    0.06416 60.00727   6.394 0.0000000268 ***
scale(sumsev_threat_t1)                      -0.14377    0.08448 59.60959  -1.702       0.0940 .  
scale(test_av_t1t2)                           0.45701    0.09573 65.30837   4.774 0.0000105802 ***
scale(test_per_cent):scale(sumsev_threat_t1) -0.08779    0.06161 60.70625  -1.425       0.1593    
scale(sumsev_threat_t1):scale(test_av_t1t2)   0.26256    0.11252 68.71931   2.333       0.0226 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) sc(__1) s(__12 s(__):
scl(tst_p_) -0.019                             
scl(sms__1)  0.068 -0.023                      
scl(ts__12) -0.032 -0.025 -0.199               
s(__):(__1) -0.018  0.088 -0.024  -0.033       
s(__1):(__1 -0.209 -0.033 -0.128   0.213 -0.051

#confint.merMod(mod2_test_threat_boys, method = "boot")
# higher threat severity in boys________________________________________________
mod2_test_threat_boys_hi <-
   lmer(
    cort_clean ~
      scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(test_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_boys
  )
summary(mod2_test_threat_boys_hi)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) + scale(test_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) -  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 198.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.40026 -0.48992  0.08387  0.64318  1.62534 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.06621  0.2573  
 Residual             0.19860  0.4456  
Number of obs: 123, groups:  ELS_ID, 64

Fixed effects:
                                                                                                                   Estimate
(Intercept)                                                                                                        -1.69066
scale(test_per_cent, scale = FALSE)                                                                                 0.44216
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     -0.04776
scale(test_av_t1t2, scale = FALSE)                                                                                  0.85459
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) -0.06260
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)   0.16212
                                                                                                                   Std. Error
(Intercept)                                                                                                           0.07810
scale(test_per_cent, scale = FALSE)                                                                                   0.12720
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                        0.02806
scale(test_av_t1t2, scale = FALSE)                                                                                    0.19303
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))    0.04393
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)     0.06948
                                                                                                                         df
(Intercept)                                                                                                        60.79473
scale(test_per_cent, scale = FALSE)                                                                                62.05761
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     59.60959
scale(test_av_t1t2, scale = FALSE)                                                                                 67.98183
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) 60.70625
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  68.71931
                                                                                                                   t value
(Intercept)                                                                                                        -21.648
scale(test_per_cent, scale = FALSE)                                                                                  3.476
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                      -1.702
scale(test_av_t1t2, scale = FALSE)                                                                                   4.427
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))  -1.425
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)    2.333
                                                                                                                               Pr(>|t|)
(Intercept)                                                                                                        < 0.0000000000000002
scale(test_per_cent, scale = FALSE)                                                                                            0.000935
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                                 0.093982
scale(test_av_t1t2, scale = FALSE)                                                                                            0.0000355
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))             0.159270
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)              0.022563
                                                                                                                      
(Intercept)                                                                                                        ***
scale(test_per_cent, scale = FALSE)                                                                                ***
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     .  
scale(test_av_t1t2, scale = FALSE)                                                                                 ***
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))    
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F-sn=T s(__1s=F ss=Fs=F-sn=T
s(__,s=FALS  -0.039                                         
I((s=F-sn=T   0.737 -0.032                                  
s(__12,s=FA  -0.254 -0.062   -0.206                         
ss=Fs=F-sn=T -0.029  0.725   -0.024    -0.055               
Is=F-sn=Ts=F -0.231 -0.057   -0.128     0.818   -0.051      

mod2_test_threat_boys_hi_int <- summary(mod2_test_threat_boys_hi)$coefficients[1]
mod2_test_threat_boys_hi_slp <- summary(mod2_test_threat_boys_hi)$coefficients[2]
# lower threat severity in boys________________________________________________
mod2_test_threat_boys_lo <-
   lmer(
    cort_clean ~
      scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(test_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_boys
  )
summary(mod2_test_threat_boys_lo)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) + scale(test_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) +  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_boys

REML criterion at convergence: 198.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.40026 -0.48992  0.08387  0.64318  1.62534 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.06621  0.2573  
 Residual             0.19860  0.4456  
Number of obs: 123, groups:  ELS_ID, 64

Fixed effects:
                                                                                                                   Estimate
(Intercept)                                                                                                        -1.50694
scale(test_per_cent, scale = FALSE)                                                                                 0.68296
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     -0.04776
scale(test_av_t1t2, scale = FALSE)                                                                                  0.23094
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) -0.06260
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)   0.16212
                                                                                                                   Std. Error
(Intercept)                                                                                                           0.07294
scale(test_per_cent, scale = FALSE)                                                                                   0.11653
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                        0.02806
scale(test_av_t1t2, scale = FALSE)                                                                                    0.15592
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))    0.04393
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)     0.06948
                                                                                                                         df
(Intercept)                                                                                                        58.34553
scale(test_per_cent, scale = FALSE)                                                                                58.35685
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     59.60959
scale(test_av_t1t2, scale = FALSE)                                                                                 66.19008
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) 60.70625
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  68.71931
                                                                                                                   t value
(Intercept)                                                                                                        -20.661
scale(test_per_cent, scale = FALSE)                                                                                  5.861
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                      -1.702
scale(test_av_t1t2, scale = FALSE)                                                                                   1.481
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))  -1.425
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)    2.333
                                                                                                                               Pr(>|t|)
(Intercept)                                                                                                        < 0.0000000000000002
scale(test_per_cent, scale = FALSE)                                                                                         0.000000227
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                                   0.0940
scale(test_av_t1t2, scale = FALSE)                                                                                               0.1433
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))               0.1593
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)                0.0226
                                                                                                                      
(Intercept)                                                                                                        ***
scale(test_per_cent, scale = FALSE)                                                                                ***
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     .  
scale(test_av_t1t2, scale = FALSE)                                                                                    
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))    
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F+sn=T s(__1s=F ss=Fs=F+sn=T
s(__,s=FALS  -0.001                                         
I((s=F+sn=T  -0.691  0.000                                  
s(__12,s=FA   0.139 -0.006   -0.035                         
ss=Fs=F+sn=T  0.004 -0.659   -0.024     0.020               
Is=F+sn=Ts=F -0.057  0.012   -0.128    -0.702   -0.051      

mod2_test_threat_boys_lo_int <- summary(mod2_test_threat_boys_lo)$coefficients[1]
mod2_test_threat_boys_lo_slp <- summary(mod2_test_threat_boys_lo)$coefficients[2]

mod2_test_threat_girls <-
 lmer(
    scale(cort_clean) ~
      scale(test_per_cent) * scale(sumsev_threat_t1) +
      scale(test_av_t1t2) * scale(sumsev_threat_t1) +
      (1 | ELS_ID),
    data = cd_clean_girls
  )
summary(mod2_test_threat_girls)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_per_cent) * scale(sumsev_threat_t1) +  
    scale(test_av_t1t2) * scale(sumsev_threat_t1) + (1 | ELS_ID)
   Data: cd_clean_girls

REML criterion at convergence: 367.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.51504 -0.64320 -0.00077  0.52372  1.88116 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1859   0.4311  
 Residual             0.6722   0.8199  
Number of obs: 133, groups:  ELS_ID, 73

Fixed effects:
                                              Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)                                   0.008466   0.088423 66.515888   0.096 0.924014    
scale(test_per_cent)                          0.180267   0.072494 66.239074   2.487 0.015423 *  
scale(sumsev_threat_t1)                       0.034255   0.093726 65.609394   0.365 0.715929    
scale(test_av_t1t2)                           0.369763   0.102718 68.159150   3.600 0.000599 ***
scale(test_per_cent):scale(sumsev_threat_t1) -0.175229   0.075262 62.992713  -2.328 0.023123 *  
scale(sumsev_threat_t1):scale(test_av_t1t2)  -0.066352   0.076732 64.356233  -0.865 0.390400    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) sc(__) sc(__1) s(__12 s(__):
scl(tst_p_) -0.030                             
scl(sms__1)  0.059  0.021                      
scl(ts__12) -0.049 -0.003  0.046               
s(__):(__1)  0.012 -0.066 -0.021   0.036       
s(__1):(__1  0.120  0.024  0.252  -0.320 -0.030

#confint.merMod(mod2_test_threat_girls, method = "boot")
# higher threat severity in girls________________________________________________
mod2_test_threat_girls_hi <-
   lmer(
    cort_clean ~
      scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(test_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_girls
  )
summary(mod2_test_threat_girls_hi)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) + scale(test_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) -  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_girls

REML criterion at convergence: 241.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.51504 -0.64320 -0.00077  0.52372  1.88116 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.07199  0.2683  
 Residual             0.26039  0.5103  
Number of obs: 133, groups:  ELS_ID, 73

Fixed effects:
                                                                                                                   Estimate
(Intercept)                                                                                                        -1.45668
scale(test_per_cent, scale = FALSE)                                                                                 0.01254
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                      0.01012
scale(test_av_t1t2, scale = FALSE)                                                                                  0.52159
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) -0.20716
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  -0.05416
                                                                                                                   Std. Error
(Intercept)                                                                                                           0.08252
scale(test_per_cent, scale = FALSE)                                                                                   0.25140
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                        0.02770
scale(test_av_t1t2, scale = FALSE)                                                                                    0.18347
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))    0.08898
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)     0.06264
                                                                                                                         df
(Intercept)                                                                                                        65.68829
scale(test_per_cent, scale = FALSE)                                                                                61.25590
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                     65.60939
scale(test_av_t1t2, scale = FALSE)                                                                                 63.76436
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) 62.99271
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  64.35623
                                                                                                                   t value
(Intercept)                                                                                                        -17.653
scale(test_per_cent, scale = FALSE)                                                                                  0.050
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                       0.365
scale(test_av_t1t2, scale = FALSE)                                                                                   2.843
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))  -2.328
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)   -0.865
                                                                                                                              Pr(>|t|)
(Intercept)                                                                                                        <0.0000000000000002
scale(test_per_cent, scale = FALSE)                                                                                             0.9604
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                                  0.7159
scale(test_av_t1t2, scale = FALSE)                                                                                              0.0060
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))              0.0231
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)               0.3904
                                                                                                                      
(Intercept)                                                                                                        ***
scale(test_per_cent, scale = FALSE)                                                                                   
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE))                                        
scale(test_av_t1t2, scale = FALSE)                                                                                 ** 
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)) *  
I(scale(sumsev_threat_t1, scale = FALSE) - sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F-sn=T s(__1s=F ss=Fs=F-sn=T
s(__,s=FALS  -0.009                                         
I((s=F-sn=T   0.746 -0.001                                  
s(__12,s=FA   0.186  0.020    0.226                         
ss=Fs=F-sn=T -0.007  0.698   -0.021     0.013               
Is=F-sn=Ts=F  0.258 -0.005    0.252     0.411   -0.030      

mod2_test_threat_girls_hi_int <- summary(mod2_test_threat_girls_hi)$coefficients[1]
mod2_test_threat_girls_hi_slp <- summary(mod2_test_threat_girls_hi)$coefficients[2]
# lower threat severity in girls________________________________________________
mod2_test_threat_girls_lo <-
   lmer(
    cort_clean ~
      scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      scale(test_av_t1t2, scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) +
      (1 | ELS_ID),
    data = cd_clean_girls
  )
summary(mod2_test_threat_girls_lo)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cort_clean ~ scale(test_per_cent, scale = FALSE) * I(scale(sumsev_threat_t1,  
    scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) + scale(test_av_t1t2,  
    scale = FALSE) * I(scale(sumsev_threat_t1, scale = FALSE) +  
    sd(sumsev_threat_t1, na.rm = TRUE)) + (1 | ELS_ID)
   Data: cd_clean_girls

REML criterion at convergence: 241.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.51504 -0.64320 -0.00077  0.52372  1.88116 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.07199  0.2683  
 Residual             0.26039  0.5103  
Number of obs: 133, groups:  ELS_ID, 73

Fixed effects:
                                                                                                                   Estimate
(Intercept)                                                                                                        -1.49931
scale(test_per_cent, scale = FALSE)                                                                                 0.88506
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                      0.01012
scale(test_av_t1t2, scale = FALSE)                                                                                  0.74972
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) -0.20716
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  -0.05416
                                                                                                                   Std. Error
(Intercept)                                                                                                           0.07780
scale(test_per_cent, scale = FALSE)                                                                                   0.26864
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                        0.02770
scale(test_av_t1t2, scale = FALSE)                                                                                    0.25199
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))    0.08898
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)     0.06264
                                                                                                                         df
(Intercept)                                                                                                        66.42664
scale(test_per_cent, scale = FALSE)                                                                                67.54860
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                     65.60939
scale(test_av_t1t2, scale = FALSE)                                                                                 68.42245
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) 62.99271
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)  64.35623
                                                                                                                   t value
(Intercept)                                                                                                        -19.271
scale(test_per_cent, scale = FALSE)                                                                                  3.295
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                       0.365
scale(test_av_t1t2, scale = FALSE)                                                                                   2.975
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))  -2.328
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)   -0.865
                                                                                                                               Pr(>|t|)
(Intercept)                                                                                                        < 0.0000000000000002
scale(test_per_cent, scale = FALSE)                                                                                             0.00157
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                                  0.71593
scale(test_av_t1t2, scale = FALSE)                                                                                              0.00404
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))              0.02312
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)               0.39040
                                                                                                                      
(Intercept)                                                                                                        ***
scale(test_per_cent, scale = FALSE)                                                                                ** 
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE))                                        
scale(test_av_t1t2, scale = FALSE)                                                                                 ** 
scale(test_per_cent, scale = FALSE):I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)) *  
I(scale(sumsev_threat_t1, scale = FALSE) + sd(sumsev_threat_t1, na.rm = TRUE)):scale(test_av_t1t2, scale = FALSE)     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) s(__,s=F Is=F+sn=T s(__1s=F ss=Fs=F+sn=T
s(__,s=FALS  -0.042                                         
I((s=F+sn=T  -0.708  0.029                                  
s(__12,s=FA   0.006 -0.038   -0.100                         
ss=Fs=F+sn=T  0.025 -0.742   -0.021     0.041               
Is=F+sn=Ts=F -0.104  0.037    0.252    -0.748   -0.030      

mod2_test_threat_girls_lo_int <- summary(mod2_test_threat_girls_lo)$coefficients[1]
mod2_test_threat_girls_lo_slp <- summary(mod2_test_threat_girls_lo)$coefficients[2]

plot_mod2_test <-
  tibble(
    Sex = c("Girls", "Boys", "Girls", "Boys"),
    Threat = c("Higher (+1SD)", "Higher (+1SD)", "Lower (-1SD)", "Lower (-1SD)"),
    Intercept = c(
      mod2_test_threat_girls_hi_int, 
      mod2_test_threat_boys_hi_int, 
      mod2_test_threat_girls_lo_int, 
      mod2_test_threat_boys_lo_int
    ),
    Slope = c(
      mod2_test_threat_girls_hi_slp, 
      mod2_test_threat_boys_hi_slp, 
      mod2_test_threat_girls_lo_slp, 
      mod2_test_threat_boys_lo_slp
    )
  )
# produce a plot for legend only_______________________________________________
plot_mod2_test_threat_legend <- 
  cd_clean_boys %>% 
  ggplot(aes(test_per_cent, cort_clean)) +
  geom_point(alpha = 1/2) +
  geom_abline(
    data = plot_mod2_test %>% filter(Sex == "Boys"),
    aes(
      intercept = Intercept,
      slope = Slope,
      color = Threat
    ),
    size = 2
  ) +
  scale_x_continuous(breaks = seq.int(-4, 4, 1)) +
  scale_color_manual(
    values = c("darkred", "darkblue")
  ) +
  expand_limits(y = c(-1, -4)) +
  expand_limits(x = c(-1, 1)) +
  theme_apa(base_size = 20, box = TRUE) +
  theme(
    plot.title = element_text(hjust = .5)
  ) +
  labs(
    title = "Boys",
    shape = NULL,
    size = NULL,
    color = "Threat severity",
    x = NULL,
    y = "log Cortisol (µg/dL)"
  )  
# plot for boys________________________________________________________________
plot_mod2_test_threat_boys <- 
  cd_clean_boys %>% 
  ggplot(aes(test_per_cent, cort_clean)) +
  geom_point(alpha = 1/2, aes(size = sumsev_threat_t1, show.legend = FALSE)) +
  geom_abline(
    data = plot_mod2_test %>% filter(Sex == "Boys"),
    aes(
      intercept = Intercept,
      slope = Slope,
      color = Threat
    ),
    size = 2
  ) +
  scale_x_continuous(breaks = seq.int(-4, 4, 1)) +
  scale_color_manual(
    values = c("darkred", "darkblue")
  ) +
  expand_limits(y = c(-1, -4)) +
  expand_limits(x = c(-1, 1)) +
  theme_apa(base_size = 20, box = TRUE) +
  theme(
    plot.title = element_text(hjust = .5)
  ) +
  labs(
    title = "Boys",
    shape = NULL,
    size = NULL,
    color = "Threat severity",
    x = NULL,
    y = NULL
  )  

Ignoring unknown aesthetics: show.legend

ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/test_cort_threat_boys.png",
  height = 6, width = 10
)
# plot for girls________________________________________________________________
plot_mod2_test_threat_girls <- 
  cd_clean_girls %>% 
  ggplot(aes(test_per_cent, cort_clean)) +
  geom_point(alpha = 1/2, aes(size = sumsev_threat_t1), show.legend = FALSE) +
  geom_abline(
    data = plot_mod2_test %>% filter(Sex == "Girls"),
    aes(
      intercept = Intercept,
      slope = Slope,
      color = Threat
    ),
    size = 2
  ) +
  scale_x_continuous(breaks = seq.int(-4, 4, 1)) +
  scale_color_manual(
    values = c("darkred", "darkblue")
  ) +
  expand_limits(y = c(-1, -4)) +
  expand_limits(x = c(-1, 1)) +
  theme_apa(base_size = 24, box = TRUE) +
  theme(
    plot.title = element_text(hjust = .5)
  ) +
  labs(
    title = "Girls",
    shape = NULL,
    size = NULL,
    color = "Threat severity",
    x = NULL,
    y = NULL
  )  
ggsave(
  "~/Box/lucy_king_files/ELS/cort_dhea/test_cort_threat_girls.png",
  height = 6, width = 10
)
# save together________________________________________________________________
plot_mod2_test_threat_sex <- 
  plot_grid(
    plot_mod2_test_threat_girls + theme(legend.position = "none"),
    plot_mod2_test_threat_boys + theme(legend.position = "none"), 
    align = "hv", 
    axis = "bt", 
    rel_widths = c(1, 1)
  ) 
#create common x and y labels___________________________________________________
y.grob_cort <- textGrob(
  "log Cortisol (µg/dL)",
  gp = gpar(col = "black", fontsize = 22), 
  rot = 90
)
x.grob_test <- textGrob(
  expression("Panel B: "*Delta*" log testosterone (pg/mL) (person-mean-centered)"),
  gp = gpar(col = "black", fontsize = 22)
)
#add common labels to plot______________________________________________________
plot_mod2_test_threat_sex <-
  grid.arrange(
  arrangeGrob(
    plot_mod2_test_threat_sex, 
    left = y.grob_cort, 
    bottom = x.grob_test
  )
)

# extract legend________________________________________________________________
threat_legend <- get_legend(
  # create some space to the left of the legend
  plot_mod2_test_threat_legend + 
    theme(
      legend.box.margin = margin(0, 0, 0, 12)
    )
)

Removed 25 rows containing missing values (geom_point).

# add legend back_______________________________________________________________
# the width of one plot (via rel_widths).
plot_grid(plot_mod2_test_threat_sex, threat_legend, rel_widths = c(1, .2))
# save plot_____________________________________________________________________
ggsave("~/Box/lucy_king_files/ELS/cort_dhea/test_cort_threat_sex.png", width = 12, height = 6)

# save DHEA and testosterone plots together
plot_mod2 <-
  grid.arrange(plot_mod2_threat_sex, plot_mod2_test_threat_sex)

plot_grid(
  plot_mod2,
  threat_legend,
  rows = 2,
  rel_heights = c(5, 1)
)

Argument 'rows' is deprecated. Use 'nrow' instead.

ggsave("~/Box/lucy_king_files/ELS/cort_dhea/hormones_cort_threat_sex.png", width = 12, height = 12)

Supplementary analyses

Distinguishing testosterone and DHEA

summary(mod2_dhea_test)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(cort_clean) ~ scale(test_per_cent) + scale(test_av_t1t2) +  
    scale(dhea_per_cent) + scale(dhea_av_t1t2) + (1 | ELS_ID)
   Data: cd_clean

REML criterion at convergence: 654.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.38135 -0.56043  0.04254  0.58635  2.42319 

Random effects:
 Groups   Name        Variance Std.Dev.
 ELS_ID   (Intercept) 0.1912   0.4373  
 Residual             0.5701   0.7550  
Number of obs: 253, groups:  ELS_ID, 135

Fixed effects:
                       Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)           -0.002022   0.060866 129.457790  -0.033 0.973554    
scale(test_per_cent)   0.219644   0.059762 120.646344   3.675 0.000356 ***
scale(test_av_t1t2)    0.105807   0.081793 126.573617   1.294 0.198162    
scale(dhea_per_cent)   0.176903   0.059659 121.749561   2.965 0.003641 ** 
scale(dhea_av_t1t2)    0.306472   0.082404 138.370495   3.719 0.000290 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) scl(t__) scl(t__12) scl(d__)
scl(tst_p_) -0.013                             
scl(ts__12) -0.005 -0.012                      
scl(dh_pr_) -0.007 -0.587   -0.007             
scl(dh__12)  0.001 -0.004   -0.577      0.022  

Effects of non-threat-related ELS

write_csv(cd_clean, "~/Box/lucy_king_files/ELS/cort_dhea/data_final/els_cort_hormones_final.csv")

cannot open file '/Users/lucyking/Box/lucy_king_files/ELS/cort_dhea/data_final/els_cort_hormones_final.csv': No such file or directoryError in open.connection(path, "wb") : cannot open the connection